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ABSTRACT 

\ 

We  study  finite  difference  approximations  to  weak  solutions  of  the  Cauchy 
problem  for  hyperbolic  systems  of  conservation  laws  in  one  space  dimension. 

We  establish  stability  in  the  total  variation  norm  and  convergence  for  a  class 
of  hybridized  schemes  which  employ  the  random  choice  scheme  together  with 
perturbations  of  classical  conservative  schemes.  We  also  establish  partial 
stability  results  for  classical  conservative  schemes.  Our  approach  is  based 
on  an  analysis  of  finite  difference  operators  on  local  and  global  wave 
conf igurations . 
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SIGNIFICANCE  AND  EXPLANATION 


He  are  ooncerned  with  the  numerical  computation  of  shock  waves  using 
finite  difference  schemes.  Specifically,  we  study  problems  concerning  the 
stability  and  convergence  of  finite  difference  approximations  and  problems  of 
describing  the  propagation  of  physical  and  numerical  waves . 
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FINITE  DIFFERENCE  SCHEMES  FOR  CONSERVATION  LAWS 


Ronald  J.  DiPerna* 


1.  INTRODUCTION. 

We  are  concerned  with  finite  difference  approximations  to  weak  solutions  of 
the  Cauchy  problem  for  hyperbolic  systems  of  conservation  laws  of  the  form 

(1.1)  u  +  f(u)  =  0,-“<x<°°. 

t  x 

Here  the  solution  u  =  u(x,t)  takes  on  values  in  Rn  and  f  is  a  smooth 
nonlinear  mapping  from  Rn  to  Rn.  We  assume  that  the  system  is  strictly 
hyperbolic  in  the  sense  that  the  Jacobian  matrix  7f(u)  has  n  real  and 
distinct  eigenvalues 

X(u)  <  X_(u)  <  ...  <  X  (u)  , 

12  n 

and  we  require  that  each  eigenvalue  X^  is  either  genuinely  nonlinear  or 
linearly  degenerate  in  the  sense  of  Lax  [20]  ,  i.e.  either 

(1.2)  r.*V.^0  or  r.'VX.HO 

3  3  3  3 

for  each  index  j  where  r^  =  r^(u)  denotes  the  corresponding  right 
eigenvector  of  Vf(u).  Systems  with  this  structure  arise  in  several  branches  of 
continuum  mechanics:  fluid  dynamics,  MHD,  elasticity,  etc. 

Experience  with  (1.1)  has  indicated  that  the  space  BV  of  functions  of 
bounded  variation  provides  a  natural  setting  for  the  solution  operator.  It  is 
known  for  example  that  if  the  initial  data  Uq(x)  lie  in  a  small  neighborhood 
of  a  fixed  state  u  e  Rn  and  have  small  total  variation,  then  a  subsequence  of 
the  family  of  difference  approximations  u(x,t,Ax)  generated  by  the  random 
choice  method  of  Glimm  converges  pointwise  a.e.  to  a  solution  u  [12]. 
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Moreover,  the  entire  family  is  stable  in  the  total  variation  norm  in  the  sense 
that 

(1.3)  TVu(*,t,Ax)  <  const  TVUg(*)  , 

where  the  constant  depends  only  on  f  and  u;  corresponding  estimates  hold  in 
the  limit  for  u.  In  the  case  of  data  with  large  total  variation,  analogous 
stability  and  convergence  results  for  the  random  choice  method  have  been 
obtained  for  certain  special  systems  [1,  7,  8,  17,  24,  25,  28,  29].  These 
results  motivate  the  problem  of  determining  the  mechanisms  which  induce  or 
preclude  stability  in  the  total  variation  norm  for  standard  finite  difference 
schemes,  i.e.  schemes  which  are  conservative  in  the  Lax-Wendroff  sense  [22]. 

We  note  that  the  problem  of  establishing  stability  for  conservative 

oo  2 

difference  schemes  in  any  of  the  natural  spaces  for  (1.1),  e.g.  BV,  L  ,  L  and 
convergence  remains  open  except  in  the  case  of  first  order  accurate  methods 
applied  to  scalar  conservation  laws  [5,  30,  34]:  in  the  case  n  =  1  the 
structure  of  the  equation  induces  maximum  principles  for  the  corresponding  exact 
solution  operator  in  BV  and  iP ,  1  <  p  <  <*>;  these  maximum  principles  are 
preserved  by  the  difference  operators  of  those  schemes  which  are  precisely  first 
order  accurate.  With  regard  to  L  -stability  and  its  relationship  to  proper 
entropy  production,  we  refer  the  reader  to  the  work  of  Majda  and  Osher  (27]  on 
second  order  accurate  schemes  applied  to  scalar  equations.  In  connection  with 

OO 

the  related  role  of  L  -stability  and  action  of  the  exact  solution  operator  in 
the  weak  topology,  we  reference  the  work  of  Tartar  [32]  on  the  theory  of 
compensated  compactness,  which  contains  several  convergence  results  for  exact 
solutions  to  general  scalar  conservation  laws  and  their  associated  parabolic 
regularizations.  In  the  computational  setting,  we  refer  the  reader  to  the  work 
of  Chorln  [2,  3,  4]  on  the  implementation  of  the  random  choice  method  for 


reacting  (and  non-reacting)  gas  flow,  to  the  work  of  Glimm,  Marchesin  and 


McBryan  [14,  15,  16]  on  hybridized  approaches  involving  the  random  choice  method 


together  with  procedures  of  subgrid  refinement  and  wave  tracking,  to  Crandall 
and  Majda  [35]  on  fractional  step  methods  and  to  Engquist  and  Osher  [36]  on  one¬ 
sided  difference  approximations. 

In  this  paper  we  are  primarily  concerned  with  theoretical  aspects  of 
stability  in  the  total  variation  norm  and  convergence  for  general  finite 
difference  schemes  for  systems  of  equations.  We  note  that  stability  of  the  form 
(1.3)  for  a  family  of  approximate  solutions  guarantees  the  existence  of  a 
subsequence  converging  pointwise  a.e.,  nee  equations  of  the  form  (1.1)  link 
the  temporal  and  spatial  variations  of  u;  convergence  of  the  entire  family 
follows  from  uniqueness,  in  those  particular  circumstances  where  uniqueness  is 
available  [9].  One  may,  of  course,  entertain  growth  estimates  on  the  total 
variation  norm  which  are  uniform  in  the  mesh  length. 

We  begin  in  Section  2  by  formulating  a  new  class  K  of  difference  schemes 
which  are  conservative  in  the  Lax-Wendrof f  sense;  the  class  K  arises  from  a 
discrete  approximation  to  the  contour  integral  form  of  system  (1.1)  taken  with 
respect  to  parallelograms  having  space-like  sides  in  the  x-t  plane.  The 
standard  conservative  schemes  with  a  three-point  domain  of  dependence  can  be 
subsumed  by  K  after  introducing  a  fractional  step.  In  Section  8  we  introduce 
a  new  class  of  hybridized  schemes  which  employ  the  random  choice  method  to 
approximate  shock  waves  and  perturbations  of  a  certain  class  N  of  first  order 
accurate  schemes  in  K  to  approximate  the  continuous  regions  of  flow;  these 
hybridized  schemes  are  based  on  the  tracking  of  waves  whose  magnitudes  lie 
between  two  specified  thresholds  depending  on  the  mesh  length.  In  the  case  of 
initial  data  with  small  total  variation,  we  establish  stability  in  the  total 
variation  norm  and  pointwise  a.e.  convergence  of  the  difference  approximat ions 
generated  by  the  hybridized  schemes  applied  to  a  class  of  systems  of  two 
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equations,  cf.  Sections  10  and  11.  This  class  includes  systems  of  the  form 

(1.4)  v  +  p(w)  =0,  w  +  q{v)  =  0 

t  x  t  x 

where  p'q'  >  0,  e.g.  the  isentropic  equations  of  gas  dynamics  and  the 
equations  for  thin  elastic  beams  in  Lagrangian  coordinates.  In  Section  12,  we 
verify  that  the  solutions  constructed  by  these  hybridized  methods  satisfy  the 
entropy  condition  of  Lax  [21].  For  general  systems  of  n  equations  we  obtain 
certain  partial  results  concerning  stability  in  the  total  variation  norm  for  the 
aforementioned  subclass  N  of  first  order  accurate  conservative  schemes;  the 
subclass  N  includes  the  Lax-Friedrichs  scheme,  cf.  Section  4.  For  both  the 
conservative  and  hybridized  schemes,  the  total  variation  estimates  are  obtained 
with  the  aid  of  non-monotone  functionals  which  are  equivalent  to  the  total 
variation  norm,  cf.  Sections  4,  5  and  9. 

The  form  of  these  functionals  for  schemes  in  class  N  is  motivated  by  an 
analysis  of  the  corresponding  difference  operators  on  discrete  wave 
interactions.  In  Sections  3  and  4  we  describe  a  general  approach  to  the  problem 
of  analyzing  difference  operators  on  local  and  global  wave  conf igurations  and 
apply  it  to  the  subclass  N.  For  the  purpose  of  analyzing  the  local  action,  we 
formulate  a  working  notion  of  local  discrete  wave  interaction  which  is  based  on 
a  process  of  interpolating  elementary  waves  between  adjacent  mesh  points.  We 
then  study  the  relations  which  govern  the  magnitudes  of  the  incoming  and 
outgoing  waves  in  a  local  interaction,  cf.  Sections  3,  5  and  7.  In  order  to 
describe  the  global  action,  we  introduce  a  state  space  fi  of  global  wave 
configurations  X  and  associate,  with  a  given  scheme,  a  mapping 

M  :  J)  Jl 

such  that  each  of  the  difference  approximations  u  generated  by  the  scheme 
corresponds  to  a  discrete  trajectory  of  the  form 

{MkX0  :  k  =  0,1,2,...} 


when  Xq  represents  the  initial  data,  cf.  Section  4.  The  generic  state  X 
in  SI  consists  of  a  sequence  of  local  states  X  which  record  the  structure  of 

the  local  waves  in  u.  Specific  local  states  X  are  associated  with  a  given 

finite  difference  approximation  u  by  applying  an  interpolation  map 

I  “  R°  X  R*'  +  R^ 

to  pairs  of  values  of  u  at  adjacent  mesh  points;  the  map  I  transforms  a  pair 

(v,w)  in  Rn  x  Rn  into  a  set 

{e  (v,w),...,e  ( v , w ) } 
l  n 

of  elementary  wave  magnitudes  e^(v,w)  associated  with  the  classical  solution 
of  the  Riemann  problem  for  (1.1)  with  Riemann  data  (v,w).  For  concreteness,  we 
restrict  our  attention  in  this  paper  to  three-level  schemes  which  employ  a 
stencil  based  on  four  mesh  points;  we  note  that  the  Lax-Friedrichs  scheme  and 
the  random  choice  scheme  can  be  regarded  as  having  this  form,  while  the  Lax- 
Wendroff  scheme  can  be  viewed  as  a  composition  of  two  such  schemes.  For 
stencils  with  that  geometry,  the  associated  space  fl  consists  of  states 

V  V 

X  *  X  where  X  lies  in  R  and  assumes  the  form 

Xk  =  (<5,y,s)  . 

Here  s  represents  the  value  of  the  difference  approximation  u  at  a  specified 

mesh  point  (x,t)  while  6  =  (6  ,...,5  )  and  y  =  (Y«,...,Y  )  denote  sets  of 

in  in 

interpolated  wave  magnitudes  associated  with  the  mesh  points  immediately  to  the 
left  and  right  of  (x,t). 

In  the  framework  of  the  space  n,  we  formulate  working  definitions  of 
approximate  simple  wave  and  weakly  interacting  state  and  we  identify  several 
classes  A  of  such  states  which  are  attracting  from  the  marching  map  M  of 
schemes  in  class  N  in  the  sense  that  the  restriction  of  the  orbit  M*X  to 
Ac  runs  down-hill  with  respect  to  associated  "coercive"  potential  functionals 
Qft.  The  reduction  which  these  potentials  pR  experience  in  one  time-step 


i 


acting  on  states  X  in  Ac  exceeds  the  corresponding  increment  to  the  total 
variation  norm  due  to  numerical  and/or  physical  wave  amplification.  This 
property  yields  a  uniform  estimate  on  the  total  variation  norm  for  those 
segments  of  the  orbit  which  lie  in  Ac  and  have  an  initial  p>oint  with  small 
total  variation.  This  partial  stability  result  for  schemes  in  class  N  serves 
as  the  starting  point  of  our  analysis  of  the  hybridized  schemes. 

Preliminary  to  the  construction  of  the  potential  functionals,  we  compare  in 
Section  3  the  behavior  of  the  exact  solution  with  that  of  the  difference 
operators  in  class  N  on  local  wave  interactions.  In  the  setting  of  £1,  a 
discrete  interaction  consists  of  a  collection  of  pairs  of  mesh  points  together 
with  the  corresponding  sets  of  interpolated  wave  magnitudes.  In  the  case  of 
waves  with  small  maanitude,  a  Taylor  expansion  of  the  equations  relating  the 
incoming  and  outgoing  magnitudes  produces  a  set  of  dominant  terms  with  a  fairly 
clear  numerical  interpretation;  a  comparison  with  the  corresponding  expansion 
derived  by  Glimm  [12]  for  random  choice  interactions  or,  what  is  the  same  up  to 
quadratic  terms,  exact  wave  interactions  reveals  several  numerical  mechanisms 
which  are  absent  in  the  exact  solution  operator.  As  a  preliminary  step  in  the 
direction  of  classifying  the  numerical  modes  of  wave  propagation,  we  discuss  in 
the  setting  of  class  N  several  numerical  processes  which  we  refer  to  as  self¬ 
interaction,  splitting  and  incomplete  cancellation,  cf.  Section  3.  These 
processes  are  reflected  in  the  structure  of  the  potential  functionals  Qft. 

The  motivation  for  a  general  study  of  potential  functionals  in  the  context 
of  conservative  difference  schemes  is  the  following.  We  recall  that,  for  the 
exact  solution  operator,  wave  interactions  typically  increase  wave  magnitudes: 
an  exact  solution  u(x,t)  to  a  system  of  equations  generally  admits  a  countable 
set  of  times  tn  such  that 

lim  TVu ( • , t )  >  TVu( • ,t  )  . 

tit  n 

n 
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On  the  other  hand,  in  the  setting  of  the  random  choice  method,  Glimm 
demonstrated  that  a  potential  for  wave  interaction  can  be  attributed  t<~  _ach 
wave  configuration  through  a  quadratic  functional  Qrc(u)  in  such  a  way  that 
all  weak  interactions  reduce  Qrc  more  than  they  augment  the  total  variation 
norm  [12].  The  potential  Qrc  is  quadratic  in  the  sense  that 

c. (TVu)2  <  Q  (u)  <  c_(TVu)2  , 

1  rc  2 

while  non-increasing  and  compensating  in  the  sense  that 

Q  {u(*,t,Ax)}  and  F(u)  =  TVu(*,t,Ax)  +  Q  {u(*,t,Ax)} 
rc  rc 

are  both  nonincreasing  functions  of  time,  if  the  initial  data  of  the  random 
choice  approximations  u(x,t,Ax)  have  small  total  variation.  The  structure 
of  Qrc  is  discussed  in  Section  7.  The  stability  estimate  (1.3)  follows  from 
the  equivalence  of  F  and  TV  on  small  data. 

These  results  motivate  the  problem  of  constructing  potentials  for  standard 
difference  schemes  which  compensate  for  both  the  physical  and  numerical 
amplification  waves.  Now,  in  the  setting  of  conservative  difference  schemes  two 
new  features  arise.  The  first  is  associated  with  the  existence  of  numerical 
modes  of  wave  propagation;  it  is  not  difficult  to  show  for  example  that,  as  a 
consequence  of  augmented  wave  amplification,  there  exist  no  compensating 
potentials  which  are  monotone  and  depend,  as  Qrc  does,  only  on  the  magnitudes 
of  waves  in  a  given  configuration.  The  second  is  associated  with  the  existence 
of  a  substantial  class  of  states  X  representing  shock  profiles  which  are 
reproduced  by  the  scheme  after  a  finite  number  of  time  steps  module  a  spatial 
translation,  i.e. 

M^X  =  X  ;  X  =  {Xk+q}  . 
q  a 

Clearly  any  translation  invariant  monotone  function  must  be  constant  along  the 
entire  orbit  corresponding  to  each  shock  profile  X.  We  note  that  the  existence 
of  shock  profiles  for  a  broad  class  of  conservative  schemes  has  been  established 


by  Majda  and  Ralston  [28]  in  the  context  of  systems  of  equations  and  by  Jennings 


[18]  in  the  context  of  scalar  equations;  numerical  evidence  has  indicated  that 
such  states  are  stable. 

The  existence  of  numerical  wave  amplification  and  shock  profiles  motivates 
the  study  of  functionals  which  appeal  to  the  geometric  strucuture  of  wave 
configurations  in  addition  to  information  on  their  individual  wave  magnitudes 
and  which  are  non-monotone  when  restricted  to  orbits  MX.  In  particular  it 
leads  one  to  ask  if  there  exist  special  classes  A  of  states  containing  the 
shock  profiles  and  corresponding  potentials  which  are  coercive  on  Ac  in 

the  sense  that 

(1.5)  Qft(MX)  -  Qft(X)  <  -Aa(X) 

Q 

if  X  e  A  ,  where  A  (X)  denotes  the  distance  from  X  to  A  in  some  metric 

A 

on  (1  and  compensating  on  Ac  in  the  sense  that 

(1.6)  TV (MX )  +  Q  (MX )  <  TV(X )  +  Q  (X) 

A  A 

Q 

if  X  e  A  .  If  such  potentials  exist  and  if  the  scheme  under  consideration  is 

in  fact  convergent,  one  might  expect  that  the  structure  of  states  in  A  and/or 

the  coercive  behavior  of  QR  would  permit  only  a  mild  growth  independent  of  the 

mesh  length  for  the  functional  F  =  TV  +  along  the  entire  orbit. 

In  Section  5  we  construct  potentials  for  the  class  N  where  the  role 

of  A  is  played  by  certain  classes  of  approximate  simple  waves  and  by  certain 

classes  of  weakly  interacting  states.  The  functionals  are  not  monotone 

when  restricted  to  the  orbit  M  X  but  do  exhibit  a  rather  strong  coercive 

behavior  on  Ac,  satisfying  inequalities  of  the  form  (1.5)  and  (1.6)  on  Ac  ; 

2 

here  the  quantity  A  (X)  does  not  arise  exactly  as  the  square  of  a  distance 

A 

from  X  to  A  in  a  fixed  metric  on  but  rather  involves  the  square  of  a 

variable  distance  from  X  to  a  subclass  of  A.  We  conjecture  that,  along  the 
entire  orbit,  the  corresponding  functionals  F  =TV  +  for  schemes  in  class 
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N  experience  only  a  mild  growth  independent  of  the  mesh  length.  For  the 
hybridized  schemes  we  show  that  this  is  in  fact  the  case  by  appealing  to  the 
improved  resolution  of  local  wave  interactions  which  hybridization  affords. 


2.  CONSERVATIVE  DIFFERENCE  SCHEMES. 

In  this  section  we  formulate  a  new  class  of  difference  schemes  which  are 

conservative  in  the  Lax-Wendroff  sense  [22];  the  motivation  is  the  following. 

00 

Suppose  u  =  u(x,t)  is  a  distributional  solution  in  BV  n  L  to  a  system  of 
conservation  laws  (1.1):  the  vector-field  (f(u),u)  is  divergence-free  in  the 
sense  that  the  sum  of  the  measures  u^.  and  fx  vanishes  on  all  Borel  sets  R, 

(2.1)  {u  +  f(u)  }(B)  =  0  . 

t  x 

Green's  theorem  for  measures  [10,  33]  yields  an  equivalent  formulation  by 
requiring  that  the  integral  of  the  normal  component  of  (f,u)  vanish  for  all 
piecewise  smooth  closed  contours  C: 

(2.2)  /  v  u  +  v  f (u)ds  =  0 

C  * 

where  v  =  (v^.v^)  denotes,  for  concreteness,  the  outward  unit  normal  to  C 
and  ds  the  element  of  arc  length.  Indeed,  (2.2)  implies  (2.1)  provided  only 
that  C  lies  within  a  substantial  class  of  contours,  for  example, 
parallelograms  with  sides  parallel  to  two  fixed  directions. 

Classical  conservative  schemes  correspond  to  a  discrete  approximation  of 
(2.1)  with  C  taken  as  the  boundary  of  a  rectangle  with  sides  parallel  to  the 
axes:  for  example,  the  standard  conservative  schemes  with  a  three-point  domain 

of  dependence  employ  a  grid  with  mesh  points  of  the  form  (idx,jdt),  i  and 
j  arbitrary  integers,  and  generate  the  value  of  the  difference  approximation, 
say  u,  at  a  typical  mesh  point  (x,t)  in  terms  of  the  three  known  values 


immediately  below 


(2.3)  u(x,t)  =  ^>{u(x  -  Ax,t  -  At),  u( x ,t  -  At),  u(x  +  Ax,t  -  At)}  ; 

the  generating  function  ip  depends  on  the  choice  of  mesh  lengths  Ax  and  At 
and  is  derived  from  a  discrete  approximation  to  (2.1)  with  C  taken  as  the 
boundary  of  the  rectangle 

{(x,t)  :  iAx  <  x  <  (i  +  l)Ax,  jAt  <  t  <  (j  +  1 )At}  . 

In  this  section  we  shall  describe  a  class  of  conservative  schemes  based  on 
a  grid  having  a  diamond-shaped  stencil,  i.e.  mesh  points  of  the  form  (iAx, jAt) 
where  i  and  j  are  integers  such  that  i  +  j  is  even.  We  begin  by 
describing  a  class  of  three-level  schemes  where  generating  function  <j>  is 
derived  from  a  discrete  approximation  to  (2.1)  with  C  taken  as  the  boundary  of 
a  rhombus  D  with  vertices  at  three  time-levels  of  the  form 

(2.4)  n  =  (x,t),  s  *  (x,t-2At),  w  =  (x-Ax,t-At),  e  =  (x+Ax,t-At)  , 

where  (x,t)  is  a  typical  mesh  point:  the  value  uR  of  the  difference 

approximation  at  the  north  vertex  n  is  generated  from  known  values  at  the 

west,  south  and  east  vertices  by  a  formula  of  the  form 

u  =  <(>(u  ,u  ,u  )  , 
n  w  s  e 

where  4>  depends  on  the  ratio  of  mesh  lengths  Ax/At.  To  be  precise,  let 

a  =  (a  ,a  )  and  8  =  (8  ,8  )  denote  respectively  the  outward  unit  normals  to 
t  x  t  x 

the  ne  side  (northeast)  and  the  wn  side  of  the  rhombus  D  with  verticies 
(2.3).  The  normals  a  and  8  depend  only  on  the  fixed  ratio  of  mesh  lengths 
Ax/ At.  We  introduce  two  smooth  mappings 

H  :  Rn  x  R°  >  Rn  and  G  :  r”  x  Rn  -*•  Rn 

which  respectively  reduce  along  the  diagonal  a  =  b  to  the  normal  component  of 
the  vector  (f(u),u)  in  the  directions  a  and  8: 


(2.5) 


H(a,b)  =  ata  +  axf(a)  +  p(a,b)(a  -  b) 
G(a,b)  =  8fca  +  8xf(a)  +  q(a,b)(a  -  b)  . 


-10- 


Here  p  and  q  simply  denote  arbitrary  smooth  n  x  n  matrices;  appropriate 
restrictions  will  be  placed  on  p  and  q  below.  Next,  we  introduce  the 
following  formal  approximations  to  the  four  line  integrals  whose  sum  represents 
the  contour  integral  around  the  rhombus  D: 


H (u  , u  ) As  ~  /  a  u  +  a  f(u)ds?  G(u  ,u  )As  ~  /  B^u  +  B  f(u)ds 
ne  '  t  x  nw  '  t  x 

ne  wn 

-H(u  ,u  ) As  -  /  a  u  +  a  f(u)ds;  -G(u  ,u  )As  -  /  B^u  +  B  f(u)ds  , 

w  s  J  t  x  e  s  1  t  x 

ws  se 


where  As  denotes  the  length  of  the  sides  of  D.  Summing  the  formulas  above 

yields  a  formal  approximation  of  the  contour  integral  (2.2)  with  C  replacing 

by  D;  the  sum  becomes  an  equation  of  the  form 

(2.6)  H(u  , u  )  +  G(u  , u  )  -  H(u  ,u  )  -  G(u  ,u  )  =  0  , 

ne  nw  ws  es 

after  the  common  coefficient  As  is  factored  out.  Finally  we  assume  that  (2.6) 
can  be  solved  for  uR  in  terms  of  the  remaining  variables  yielding  a  smooth 
generating  function 


For  local  purposes,  solvability  is  guaranteed  by  requiring  that  the  matrix 

H  (a, a)  +  G  (a, a) 

be  invertible.  For  the  centered  schemes  described  above  we  have 

a  =  B  ,  a  =  -B  and 
t  t  x  x 

H  (a, a)  +  G  (a, a)  =  2a  I  +  p(a,a)  +  q(a,a) 

cl  8  u 

which  is  invertible  if  for  example  p  +  q  is  small  on  the  diagonal. 

In  a  similar  fashion,  one  can  construct  non-centered  conservative  schemes 
based  on  a  discrete  approximation  to  (2.2)  with  C  taken  as  the  boundary  of  a 
parallelogram  with  vertices  of  the  form 

n  =  (x,t),  s  =  (x  +  eAx,t  -  2At) ,  w  *  (x  +  Ax,t  -  At),  e  =  (x  -  6Ax,t  -  At)  . 


More  generally,  one  can  construct  a  class  K  of  multi-level  conservative 


schemes  based  on  a  discrete  approximation  to  (2.2)  with  C  taken  as  the 

boundary  of  a  parallelogram  intersecting  several  time-levels.  However,  we  shall 

restrict  our  attention  for  concreteness  to  the  subclass  of  centered  three-level 

schemes  based  a  stencil  with  vertices  of  the  form  (2.4).  We  note  that  this 

subclass  contains  several  classical  schemes.  The  leap-frog  scheme  is  formed 

from  an  arithmetic  average, 

H ( a ,b)  =  afc(a  +  b)/2  +  a  {f(a)  +  f(b)}/2 

G(a,b)  =  fl  (a  +  b)/2  +  6  {f(a)  +  f(b)}/2  . 
t  x 

The  Lax-Friedrichs  scheme  is  obtained  by  eliminating  the  dependence  on  ug , 
i.e.  by  taking  p  =  q  =  0.  The  general  three-point  conservative  scheme  (2.3) 
can  be  regarded  as  a  composition  of  two  schemes  in  this  subclass  by  introducing 
a  fractional  step  to  produce  a  nine-point  stencil  having  mesh  points  on  five 
levels  of  the  form 

(x,t) , ( x , t  ±  At) , (x  ±  Ax,t) , (x  ±  ^  Ax,t  +  j  At) , (x  t  ~  Ax, t  -  j  At)  . 

Indeed,  the  standard  two-step  Lax-Wendroff  scheme  is  already  in  this  form  since 
it  can  be  regarded  as  a  composition  of  the  Lax-Friedrich  and  leap-frog  schemes. 

3.  DISCRETE  WAVES. 

In  this  section  we  shall  describe  a  method  for  introducing  local  wave 
magnitudes  into  a  finite  difference  approximation.  We  shall  present  the  method 
in  the  setting  of  the  class  K1  of  three-level  schemes  with  a  centered  diamond 
shaped  stencil;  it  has  an  obvious  analogue  for  more  general  stencils.  We  shall 
also  discuss  the  equations  which  relate  the  incoming  and  outgoing  magnitudes  of 
a  local  interaction  and  compare  them  with  the  corresponding  equations  for  the 
exact  solution  and  the  random  choice  method. 

A  pattern  of  local  wave  magnitudes  can  be  associated  with  a  difference 


approximation  defined  on  a  grid  by  using  the  classical  solution  of  the  Riemann 


problem  [20] ,  i.e.  the  initial  value  problem  with  data  of  the  form 
uQ(x)  =  u  if  x  <  0,  uQ(x)  =  u+  if  x  >  0  , 
where  u-  and  u+  denote  elements  of  Rn.  We  recall  that  the  exact  solution 
operator  resolves  Riemann  data  (u_,u+)  into  a  similarity  solution  u  =  u(x/t) 
consisting  of  n  +  1  constant  states  u_. ,  j  =  1,2,...,n  +  1,  with  adjacent 
constant  states  separated  by  either  a  j-shock  wave  or  a  centered  j-rarefaction 
wave  [20].  Here  u1  =  u~  and  un+,j  =  u+.  In  the  standard  fashion,  we  take  the 
magnitude  of  a  j-shock  wave  separating  states  u^  and  Uj  +  1  to  be  the  negative 
of  the  distance  from  u^  to  Uj+i  along  the  j-shock  wave  curve  through 
and  the  magnitude  of  a  centered  j-rarefaction  wave  separating  states  u^  and 
Uj+.j  to  the  (positive)  distance  from  u^  to  uj+i  along  the  j-rarefaction 
wave  curve  through  u^.  For  example,  the  solution  u  to  the  Riemann  problem 

for  a  system  of  two  equations  might  consist  of  a  1 -shock  x  =  at  separating 
states  u~  and  u2  together  with  a  centered  2-rarefaction  wave  separating  u2 

and  u+,  i.e.  u  =  u”  in  <  x/t  <  a,  u  =  in  o  <  x/t  <  ^2^U2*'  in 

u  *=  u+  in  X2(u+)  <  x/t  <  00  and  the  section  X2(u2)  <  x^t  <  fon,,s  a 

centered  2-rarefaction  wave.  In  general  we  shall  denote  by 

=  e_j(u  ,u+) 

the  magnitude  of  the  j-wave  in  the  solution  of  the  Riemann  problem  with  data 
(u~,u+) . 

We  shall  restrict  our  attention  to  the  class  K1  of  three-level  schemes 
in  K  based  on  a  grid  having  a  centered  diamond-shaped  stencil  with  vertices 
of  the  form  (2.4).  Consider  a  corresponding  difference  approximation 
u  =  u(x,t,Ax,At)  with  small  oscillation  and  suppose  that  the  Courant- 
Friedrichs-Levy  condition  is  satisfied,  i.e. 

Ax/A t  >  max{|Xj(v)|  :  j  =  1,2,... ,n}  , 
where  the  maximum  is  taken  over  a  set  in  Rn  containing  the  range  of  the 


difference  approximation  u.  We  shall  associate  with  u  a  pattern  of  wave 
magnitudes  by  interpolating  the  solution  of  the  Riemann  problem  between  adjacent 
mesh  points  in  the  following  fashion:  associate  with  each  pair 

=  {(x,t)  ,<x  ±  Ax,t  ±  At)} 

of  adjacent  mesh  points,  the  set  of  magnitudes  of  the  n  elementary  waves  in 
the  solution  to  the  Riemann  problem  with  data 

u  =  utp^,  u+  =  u(P2>  . 

This  association  can  be  expressed  formally  by  the  map 

X  •  R11  x  Rn  +  Rn 

I(u  ,u+)  =  {e,(u  ,u+),...,£  (u  ,u+)}  . 
i  n 

By  a  local  wave  interaction  in  u  we  shall  mean  a  configuration  consisting  of  a 
mesh  diamond  having  verticies  n,s,e  and  w  of  the  form  (2.4)  together  with 
the  four  sets  of  wave  magnitudes  which  are  obtained  from  the  Riemann  problems 
associated  with  the  four  pairs  of  adjacent  verticies  (w,s) , ( s,e) , (w,n)  and 
(n,e)  and  which  are  denoted  as  follows: 


6  = 

, . . 

) 

=  I(u  ,u  ) 

n 

w  s 

a  — 

•  .  a  ) 

=  I(u  t u  ) 
w  n 

(a1 , . . 

n 

s  e 


n  e 


We  shall  refer  to  6_.  and  y  as  the  incoming  j-waves  and  to  a_.  and  6^  as 
the  outgoing  j-waves. 

We  note  that  one  may  also  interpolate  between  adjacent  mesh  points  even  if 
the  C-F-L  condition  is  not  satisfied.  In  this  case,  however,  it  is  not  so 
clear  how  to  interpret  the  interpolation.  On  the  other  hand,  if  the  C-F-L 
condition  holds,  it  is  meaningful  to  interpolate  between  mesh  points  on  the  same 
space-like  arc  as  well  as  between  points  on  the  same  time-level,  since  the  exact 
solution  operator  applied  to  Riemann  data  is  invariant  under  a  space-like 
rotation  of  the  x-t  plane.  Indeed,  the  C-F-L  condition  guarantees  that  the 
line  segments  ws,se,wn  and  ne  forming  the  boundary  of  a  typical  mesh  diamond 
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are  space-like  for  difference  approximates  with  small  oscillation.  We  remark 

that  one  is  free  to  imagine  a  set  of  n  elementary  waves  actually  crossing  the 

line  segment  joining  a  pair  of  adjacent  mesh  points  (p1fp2)/  waves  having  a 

location  specified  only  within  a  distance  Ip^pjl.  However,  we  emphasize  that 

by  definition  we  only  associate  a  set  of  n  wave  magnitudes  with  a  given  pair 

( p.j  ,P2 ) «  With  this  understanding,  one  can  summarize  a  discrete  wave 

interaction,  say  for  a  system  of  two  equations  such  as  Figure  1  of  Section  4 

where  the  symbols  u^ , s^ ,u^+^  and  w^  denote  the  values  of  the  difference 

approximation  at  the  corresponding  diamond  indexed  by  k. 

We  shall  begin  the  discussion  of  local  interactions  by  considering  the 

relationship  between  the  incoming  and  outgoing  waves.  It  is  clear  that  the 

outgoing  magnitudes  a  and  8  are  smooth  functions  of  the  incoming  magnitudes 

6  and  Y  and  the  local  base  state  us: 

(a, 6)  =  W(6,y,u  ) 
s 

We  note  that  6,y  and  u  uniquely  determine  u„  and  u  which  together 

s  we 

with  u_  uniquely  determine  u  ,a  and  3  through  the  smooth  generating 
"  n 

function  <j>.  A  Taylor  expansion  of  W  at  5  =  y  =  0  provides  the  dominant 
terms  in  the  laws  governing  the  interaction  of  weak  waves.  To  begin  with  we  let 

v.  =  <$,,y.)  and  a.  =  (a.,0.) 

3  13  3  3  3 

denote  the  incoming  and  outgoing  j -waves  and  write 

(3.1)  o  =  A(u  )v  +  0( | v|2) 

s 

where  A(u„)  is  a  smooth  2n  x  2n  matrix  and  v  =  (v„,...,v  ), 
s  In 

a  =  (o1 , . . . ,an) . 

For  simplicity  we  shall  restrict  our  attention  to  a  broad  subclass  of 
schemes  in  which  preserves  two  basic  properties  of  the  exact  solution 

operator.  In  this  connection,  we  first  recall  that  if  all  incoming  waves  of  an 
interaction  in  an  exact  solution  belong  to  the  same  characteristic  field,  say 
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the  jth  field,  then  the  magnitudes  of  the  outgoing  waves  of  the  kth  field. 


k  y  j  are  quadratic  with  respect  to  the  magnitudes  of  the  incoming  waves.  We 
shall,  first  of  all,  restrict  attention  to  those  schemes  in  which  preserve 

this  property,  i.e.  schemes  such  that 

(3.2)  °k  =  0^ I v j 1 2  > ,  if  \  =  0  for  k  ?  j  • 

Clearly  this  property  is  equivalent  to  the  statement  that  A(ug)  is  a 
tridiagonal  matrix  for  all  ug.  We  note  that  property  (3.2)  is  satisfied  by  the 
Lax-Friedrichs  scheme  and  the  leap-frog  scheme;  a  simple  criterion  for  (3.2)  is 
given  in  a  lemma  below  in  the  setting  of  schemes  in  and  can  easily  be 

checked  for  the  Lax-Firedrichs  and  leap-frog  schemes.  It  follows  as  a  corollary 
that  the  two-step  Lax-Wendroff  scheme  satisfies  the  natural  analogue  of  (3.2) 
since  it  can  be  viewed  as  the  composition  of  the  Lax-Friedrichs  scheme  and  the 
leap-frog  scheme;  indeed,  the  same  is  true  for  all  the  standard  variations  on 
the  Lax-Wendroff  scheme  since  they  share  the  same  linearization.  Finally,  we 
remark  that  the  random  choice  method  can  be  regarded  as  a  scheme  which  employs 
the  same  stencil  as  schemes  in  and  it  satisfies  (3.2);  the  laws  for  random 

choice  interactions  are  recalled  in  Section  7. 

Secondly,  we  recall  that  wave  interactions  in  an  exact  solutin  do  not 
augment  wave  magnitudes  by  more  than  a  quantity  which  is  quadratic  with  respect 
to  the  magnitudes  of  the  incoming  waves.  We  shall  restrict  our  attention 
further  to  those  schemes  in  which  preserve  this  property,  i.e.  schemes  such 

that 

n  n 

(3.3)  l  |a  |  +  |8.|  <  l  1 6  .  |  +  |y.|  +  0(  |vp)  . 

j-1  3  3  j=1  3  3 


The  condition  (3.3)  is  equivalent  to  the  condition  that 


(3.4) 


aj  =  (1  +  TjYj  + 

8.  =  y  ,<5  +  (1  -  t.)y.  +  0  (  |  v  | 2 )  , 

3  j  j  3  3 


-16- 


1 


where  the  coefficients  y  ,  and  T.  satisfy 

3  3 

(3.5)  0  <  y.(u  )  <  1,  0  <  T.(u  )  <  1  , 

3  s  3  s 

for  all  values  of  their  argument  ug.  Now,  it  is  not  difficult  to  show  that, 
within  the  subclass  of  schemes  in  which  satisfy  property  (3.2),  inequality 

(3.3)  holds  for  those  methods  which  are  precisely  first  order  accurate.  In  this 


connection  we  recall  that  any  scheme  which  is  consistent  with  the  equations  and 
which  has  a  smooth  generating  function  is  at  least  first  order  accurate.  Thus, 
condition  (3.3)  rules  out  second  order  accurate  methods.  In  particular,  the 
Lax-Friedrich  scheme  satisfies  (3.3)  while  the  leap-frog  and  Lax-Wendroff 
schemes  do  not.  It  is  also  simply  to  verify  that  the  random  choice  method 
satisfies  (3.3),  cf.  Section  7.  Finally,  we  remark  that  for  certain  technical 
reasons  we  shall  restrict  our  attention  to  the  subclass  N  of  schemes  in  K1 
which  satisfy  (3.2),  (3.3)  and 

(3.5)  0  <  u  (r  )  <  1,  0  <  T  ( u  )  <  1  . 

j  s  j  s 

The  lemma  below  contains  a  simple  criterion  for  membership  in  N  which  shows  in 
particular,  that  the  Lax-Friedrichs  scheme  belongs  to  N. 

Lemma.  Consider  a  scheme  in  class  K.  The  corresponding  matrix  A  is 
tridiagonal  if  and  only  if  the  matrix 

w(a)  =  p(a,a)  +  q(a,a)  , 
obtained  from  (2.3),  satisfies 

(3.6)  w(a)r^(a)  =  w_.(a)r^(a),  j  =  1,2,  ...,n  , 

where  r^  denotes  the  right  eigenvector  of  f'  associated  with  the  eigenvalue 

X..  If  (3.5)  holds  then 
3 

(3.7)  y  =  {a  +  a  X.(a)  +  w.(a)}/{o  +  a  X,(a)  +  w.(a)} 

j  t  x  3  3  t  x  3  j 

(3.8)  T.  =  +  e  A  (a)  +  w.(a)}/{c>  +  a  A. (a)  +  w.(a)} 

3  t  x  j  3  t  x  3  3 

where  a  =  a  +  8  ,  a  =  a  +  8  and  the  scheme  lies  in  N  if  and  only  if 

f  r  p  v  v  v 


the  eigenvalues  Wj(a)  are  such  that  the  quantities  specified  by  (3.7)  and 

(3.8)  lie  strictly  between  zero  and  one. 

Proof:  Consider  an  iteraction  associated  with  a  mesh  diamond  such  that  u  =  u 
-  s  e 

and  uw  lies  on  the  j-wave  curve  through  ug.  In  this  case  the  only  incoming 

wave  is  a  j-wave  crossing  the  ws  side  of  the  diamond.  Regard  ug  as  fixed 

and  uw  as  parametrized  by  arc  length  6  along  the  j-wave  curve  through  ug: 

u  =  u  (6) ,  u(0)=u 

w  w  w  s 

Substituting  u  =  u  (6)  and  u  =  u  into  (2.5),  solving  for  u  =  u  (6)  and 
ww  se  nn 

differentiating  with  respect  to  6  at  6=0  yields 

(3.9)  U  (0)  =  (H  +  G  >  1 (H  -  G  )u  (0)  , 

n  a  a  a  a  w 

where  the  coefficient  matrix  is  evaluated  at  (ug,ug).  A  simple  calculation 

shows  that  the  matrix  in  (3.9)  is  given  by 

{a  I  +  a  f'(a)  +  w(a)}  1 (a  I  +  a  f'(a)  +  w(a)} 
t  x  t  x 

where  a  =  ug.  By  considering  the  analogous  incoming  configuration  where 

=  ug  and  u0  lies  on  the  j-wave  curve  through  ug  we  obtain  a 

corresponding  equation 

u  (0)  =  (H  +  G  )_1 (G  -  H.)u  (0)  , 
n  a  a  a  b  e 

in  which  the  coefficient  matrix  is  given  by 

{o  I  +  a  f '  ( a)  +  w(a)}  1  {0  I  +  6  f'(a)  +  w(a)}  , 
t  x  t  x 

with  a  =  u  .  The  lemma  follows  from  the  fact  that 
s 

u  (0)  =  u  (0)  =  r . (u  )  . 
w  e  is 

Remarks :  It  follows  from  the  lemma  that  a  scheme  in  class  N  necessarily 

satisfies  the  C-F-L  condition.  Conversely,  if  a  scheme  in  K  has  a 
tridiagonal  matrix  A  then  the  C-F-L  implies  that  u_.  and  x  lie  between 
zero  and  one  provided  w^(a)  is  sufficiently  small.  The  latter  fact  applies  to 
the  Lax-Friedrichs  scheme  for  which  p  =  q  =  0. 


* 
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Next  we  shall  describe  several  numerical  mechanisms  of  interaction  present 


in  class  N  schemes#  For  this  purpose,  let  us  write  the  expansion  (3.2)  for  a 

class  N  scheme  in  the  form 

0.  =  A  ,  ( u  )v  +  B  .  ( u  )  ( v ,  v )  +  0<  |  v| 3 ) 
j  3  3  j  3  3 


BjfUgHv.v)  =*  l  b^v(uo)(vH/vv)  +  l  bjv(uc)(v^,v4) 


3k  s  j  k 


k=1 


kk  s  j'  j 


where  v.  =  (6  ,y  ),  c,  =  (a. ,6.)  and  b^.  are  bilinear  maps  from  R2  x  R2 
3  3  3  3  3  3  xx. 

2 

to  R  depending  on  the  local  base  state  us»  The  presence  of  numerical 
mechanisms  of  wave  interaction  is  revealed  by  the  structure  of  and  B  ^ . 

Self-interactions.  The  structure  of  the  operator  B..  shows  that  the  nonlinear 
interaction  between  characteristic  fields  in  a  difference  approximation  is 
substantially  larger  than  in  an  exact  solution.  By  way  of  example,  let  us 
consider  an  exact  solution  u  to  a  system  of  two  equations  which  consists  of 
two  interacting  weak  shocks  of  different  fields.  To  be  precise,  suppose  that  in 
a  strip  of  the  form  (0,t  ),  u  consists  of  a  2-shock  $2  and  a  1 -shock  y 
which  have  trajectories 

x  -  xQ  =  sfi ( t  -  tQ )  and  x  -  xQ  =  s^  ( t  -  tQ  ) ,  t  <  tQ 
and  which  approach  with  speeds  >  s^,;  while  in  a  strip  of  the  form 
(tg,00),  u  consists  of  a  1-shock  and  a  2-shock  0^  which  have 

trajectories 

x  -  x0  =  so(t  -  tQ)  and  x  -  xQ  =  sg(t  -  tQ ) ,  t  <  tQ 

and  which  recede  with  speeds  s  <  s„;  u  is  constant  in  each  of  the  four 

a  B 

sectors  defined  by  the  four  rays  above.  It  is  well-known  the  outgoing 
magnitudes  satisfy 


O,  =  Y,  +  0(y162);  B2  =  <52  +  0(y162)  . 

For  systems  of  n  equations,  two  interacting  weak  shocks  6.  and  y  qenerate 

J  * 


two  shocks  a.  and  3,  satisfying 
3  k 

a.  =  Y  .  +  0  ( y .  5.  ) ;  6.  =  5  +  0(y.6  ) 

3  3  3kkk  3  k 

together  with  reflected  waves 

£A  =  0(Vk> 

in  the  field,  i  *  j,  £  *  k.  Thus,  in  these  two  examples  the  wave 

magnitudes  are  conserved  in  each  characteristic  field  up  to  first  order  while 
the  amplification  of  individual  waves  and  production  of  reflected  waves  is  at 
most  on  the  order  the  product  of  the  approaching  waves.  The  same  statement  can 
be  made  for  multiple  wave  interactions.  Perhaps  the  simplest  formulation  is 
provided  by  the  laws  for  multiple  wave  interactions  in  the  randon  choice  cf. 
Section  7,  these  laws  coincide  with  those  of  the  exact  solution  up  to  and 
including  quadratic  terms.  For  the  purposes  of  the  present  discussion  we  only 


want  to  remark  that  for  the  exact  solution  operator  and  random  choice  operator 

the  diagonal  terms  of  B ^  vanish  identically.  In  constrast,  the  diagonal  terms 

of  B.  for  schemes  in  class  N,  i.e.  the  matricies  b,^(u  )  do  not  vanish  on 
3  kk  s 


any  open  set.  The  term 


(3.9) 


bkk(Us)(VV 


records  a  contribution  to  the  jth  field  from  the  self -interaction  of  waves  in 
the  kth  field.  If  k  *  j  then  the  term  (3.9)  represents  a  contribution  to  the 
production  of  a  reflected  waves  in  the  jth  field  due  to  self-interactions  in  the 
kth  field.  If  k  =  j  then  the  term  (3.9)  represents  a  contribution  to  the 


amplification  of  waves  in  the  jth  field  due  to  self-interactions  in  the  jth 


field.  For  example,  suppose  that  the  incoming  waves  of  a  discrete  interaction 

1.  v. 

belong  to  the  k  field,  i.e.  =  0,  i  *  k.  First,  we  see  that  outgoing  waves 
o_.  are  produced  on  the  jtb  field  j  *  k  satisfying 

°j  ’  i  *  *  • 
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t*  h 

Thus,  although  there  exist  no  incoming  waves  of  the  j1"0  field,  j  *  k,  there  do 
exist  outgoing  "reflected"  waves  of  the  field  arising  from  self-interactions 
in  the  k^-field.  Secondly,  the  waves  of  the  kfc^  field  are  themselves  amplified 


by  a  term  of  the  form  (3.9),  i.e. 

\  ‘  *k  '  \  *  ’k  *  bto,(V(W2  *  01 'V3'  • 

Wave  Splitting.  Consider  a  discrete  interaction  with  only  one  incoming  wave. 


say 

6  . 
3 


6^.  Up  to  linear  terms,  the  action  of  a  scheme  in  class  N  is  to  split 

into  two  waves  in  a  proportion  determined  by  the  local  base  state  ug: 

a  =  {1  -  u .  ( u  )}<5.  +  0(<52)/  0.  =  p  (u  >6.  +  0(62)  . 
j  is]  j  D  j  s  j  3 


For  a  general  interaction,  each  of  the  incoming  waves  is  split  and  then 
superimposed  up  to  linear  terms  in  a  fashion  determined  by  the  structure  of  the 
matrices  Aj.  The  process  of  wave  splitting  is  absent  in  the  exact  solution 
operator  to  systems  with  eigenvalues  of  the  form  (1.2).  We  remark,  in  passing, 
that  if  an  eigenvalue  A is  not  monotone  in  the  direction  r^  then  shocks  in 
the  exact  solution  can  be  split  spontaneously  through  interactions  with  smooth 
flow.  The  process  of  wave  splitting  is  also  absent  in  the  random  choice  method 
except  for  the  trivial  situation  where  a  rarefaction  wave  is  split  by  a  sample 
point.  In  the  random  choice  method  the  splitting  of  rarefaction  waves  is  not 
accompanied  by  any  form  of  wave  amplification. 

One  of  the  interesting  consequences  of  wave  splitting  a  conservative  scheme 
is  that  the  recession  of  waves  after  interaction  is  not  sharp.  In  the  special 
exact  solution  u  described  in  the  subsection  on  self-interactions,  the  two 
receding  shock  waves  and  from  the  boundary  of  an  identically  constant 

wake  region.  For  a  conservative  scheme  two  "receding"  shock  waves  are  split 
again  and  again  at  each  time  level,  leaving  a  wake  regions  with  waves  on  the 
same  order  as  the  primary  waves  themselves. 
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Incomplete  Wave  Cancellation.  The  well-known  persistence  of  oscillations 


observed  in  difference  approximations  generated  by  conservative  difference 
schemes  corresponds  to  an  incomplete  cancellation  experienced  by  interpolated 
waves  of  the  same  characteristic  field  but  of  different  sign,  i.e.  j-waves  with 
positive  and  negative  magnitude.  In  the  random  choice  method,  the  interaction 
of  a  j-raref action  wave  (positive  magnitude)  and  a  j-shock  wave  (negative 
magnitude)  leads  to  the  absorption  of  the  smaller  wave  by  the  larger  up  to 
linear  terms.  A  similar  statement  can  be  made  for  exact  solutions  by 
considering  the  effect  of  such  an  interaction  after  a  small  interval  of  time.  In 
contrast,  in  a  conservative  difference  scheme  the  larger  wave  only  absorbs  a 
fraction  of  the  smaller  in  typical  interactions.  An  analytical  discussion  of 
this  feature  is  postponed  until  Section  6. 

4.  GLOBAL  WAVE  CONFIGURATIONS. 

In  this  section  we  shall  describe  a  framework  for  studying  local  and  global 
interactions  in  a  finite  difference  approximation  generated  by  a  member  of  the 
class  K1  of  three-level  schemes  with  a  centered  diamond-shaped  stencil.  An 
analogous  treatment  suggests  itself  for  schemes  with  more  general  stencils.  To 
begin  with,  let  us  consider  an  arbitrary  function  u  which  is  defined  on  a  grid 
having  a  centered  diamond-shaped  stencil  with  mesh  lengths  Ax  and  At.  Fix 
two  consecutive  levels  t  =  mAt  and  t  =  (m  +  1)At  and  let  (s^)  and  (u^) 
denote  the  sequences  of  values  of  u  at  the  lower  and  upper  levels;  we  put 

s,  =  u{ ( 2k-1 ) Ax, mAt}  if  m  is  even  and  s,  =  u{2kAx,mAt}  if  m  is  odd, 
k  k 

u,  =  u{2kAx, (m+IAt)  if  m  is  even  and  u  =  { (2k+1 )Ax, (m+1 )At}  if  m  is  odd, 
k  k 

and  we  put 

=  i<vV  and  yk  =  I(sk'uk+1)  • 
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Figure  1  illustrates  the  case  n  =  2.  We  note  that  the  values  s^  together 

k  k 

with  the  interpolated  magnitudes  (6  ,  y  )  uniquely  determine  the  values  u^  on 
the  upper  level.  Thus,  complete  information  for  two  consecutive  levels  is 
carried  by  the  sequence 

X  =  {Xk},  xk  =  <«k,Yk,sk)  e  R3n  . 

v 

In  certain  circumstances  it  proves  useful  to  regard  each  local  state  X  as 

a  y. 

being  decomposed  into  local  states  of  the  jun  field  X.: 

k  -k  k 

X  —  ( 5  ,  y  . ,  s  )  . 

3  3  3k 

Now  in  the  case  where  u  arises  as  a  finite  difference  approximation  associated 

with  a  scheme  in  having  a  generating  function  <j> ,  the  process  of 

advancing  u  from  one  time-level  to  the  next  can  be  conveniently  represented  by 

the  following  map  M  defined  on  the  set  S3  of  all  such  states  X: 

M{<6k,Yk,sk)}  *  (0k  1,ak,uk) 
k  k 

where  the  outgoing  wave  magnitudes  a  and  0  are  obtained  from  <f>  and  the 
corresponding  incoming  waves  by  the  rule 

ak  =  I(Vwk),  0k  =  Kwk,uk+1) 

where 


\  =  *(VWi>  ' 

cf.  Figure  1.  Thus,  the  marching  map  M  represents  the  process  by  which 

k  k 

incoming  magnitudes  6  ,y  and  base  states  s^  determine  values  uk  on  the 
next  level,  which  in  turn  produce  values  wk  and  outgoing  magnitude  through  the 
generating  function  With  the  aid  of  such  sequences,  we  shall  identify  a 

given  finite  difference  approximation  u  with  the  discrete  trajectory  of  its 
data  XQ  under  M  in  S3: 


u 


{M* 


V 


p=0 
p=0  ' 
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Figure  1 


The  problem  is  to  prove  that  if  TVXQ  is  sufficiently  small  then  for  each 
time  T  >  0, 

TVMPXq  <  const  . 

if  p  <  T/At,  where  the  constant  is  independent  of  the  mesh  length.  Here  we 
define 

TVX  =  l  |6*|  +  |y*|  . 

j,k  3  3 

As  we  remarked  in  the  introduction,  the  strategy  is  to  study  potentials  for  wave 
interaction.  For  the  subclass  N  of  schemes  satisfying  (3.2),  (3.3)  and  (3.5), 
we  shall  construct  appropriate  potentials  in  several  steps.  The  first  is  to 
introduce  a  notion  of  approximate  simple  wave  as  follows.  Embed  a  diamond¬ 
shaped  stencil,  centered  or  not,  with  mesh  lengths  Ax  and  At  into  a  j-simple 
wave  u  =  u(x,t)  and  examine  the  relationship  between  the  wave  magnitudes 
5^  and  y^  obtained  by  applying  the  interpolation  map  I  to  the  values  of 
u  at  the  west-south  vertices  and  the  south-east  vertices  respectively: 


6  =  (6  ,...,6  )  =  I(u  ,u  )}  y  =  ( y  , . . . ,y  )  =  l(u  ,u  ) 


A  simple  calculation  reveals  a  restriction  on  the  ratio  of  the  magnitudes  of  the 


incoming  j -waves: 


6.  =  0  (u  )y ,  +  0(Ax2  +  At2)  , 

3  J  s  j 

sk  *  1k  ■  °> k  * 1  ■ 

where  0^  denotes  the  ratio  of  direction  cosines  between  the  normals 
a  and  0  associated  with  the  stencil  and  the  characteristic  ray  with  speed 


vv; 


e . (u  )  =  {a  +  a  X . (u  ) > / { 3  +  B  A,(u  )}  . 

js  t  xjs  t  xjs 


For  the  purpose  of  constructing  potential  functionals,  local  states 


X  =  (6^,y^,s)  satisfying 


6j '  V",j 


play  a  central  role.  We  shall  refer  to 

I\(s)  *  {X^  :  *  6j(s)y^} 

as  the  set  of  local  j-slmple  waves  passing  through  the  base  state  s  e  Rn  and 
we  shall  study  the  action  of  class  N  schemes  on  the  corresponding  global 
configurations 


n  s)  =  {X  :  5*  =  9.(s  W*}  . 

j  3  *  j 

One  can  regard  an  element  of  T(s)  as  consisting  of  a  chain  of  local  simple 

waves  (of  varing  index)  passing  through  a  sequence  s  =  (s^)  of  local  base 

states  s,  e  Rn. 
k 

For  schemes  in  class  N  we  shall  first  construct  a  functional  P  :  ft  ♦  R 
which  decreases  along  those  segments  of  the  orbit  M %  which  lie  in  the 

complement  of  a  neighborhood  of  the  simple  waves  in  the  sense  that 
(4.2)  P (MX )  -  P(X)  <  -cd  ?X)  +  ce  ^X) 

if  TVX  <<  1.  Here  and  below  we  shall  use  the  letter  c  to  denote  any  of 
various  positive  constants  which  depend  only  on  the  system  (1.1),  the  scheme 


under  consideration  and  the  vector  u  e  R°  in  the  neighborhood  of  which  all 

2 

analysis  takes  place.  The  quantity  d(X)  denotes  the  l  -distance  from 
k  k 

X  *  {($  iY  , s  )}  to  the  corresponding  set  of  simple  waves  Hs),  s  =  (s  ), 

K  K 

i  .e. 

<J2(X)  -  [  -  0,(s.  )yJ}2/{i  +  0.(8  )}  , 

j,k  1  3  X  k  3  K 

and  ®p(X)  denotes  the  sum  of  the  p*"*1  powers  of  the  wave  magnitudes  in  X, 

1*6* 

e  (X)  =  l  |vV  . 

P  j/k  3 

One  could,  of  course,  normalize  one  of  the  constants  c  in  (4.2)  to  equal  one. 
In  addition  P  is  quadratic  in  the  sense  that 

-c(TVX)2  <  P(X)  <  c(TVX)2  . 

These  estimates  motivate  a  study  of  a  general  class  of  neighborhoods  of  simple 
waves  of  the  form 

T  =  {X  :  d2(X)  <  me  (X)}  , 

mp  p 

which  can  be  regarded  as  consisting  of  approximate  simple  waves.  However,  for 

the  purposes  of  this  paper  the  sets  T  are  needed  only  in  the  case  where  m 

mp 

is  small  and  p  =  2  and  the  case  where  m  is  large  and  p  =  3.  The  case  p  =  3 
is  of  particular  interest  since  e3(X)  represents  the  approximate  rate  of 
entropy  production  associated  with  the  state  X. 

The  details  of  the  construction  of  P  are  postponed  until  Section  5;  we 
shall  presently  restrict  our  remarks  to  certain  qualitative  properties  of  P. 

We  begin  with  a  simple  observation  that  for  schemes  in  class  N  the  inequality 
(4.1)  implies  an  estimate  on  the  total  variation  along  those  segments  of  the 

Q 

orbit  in  T  .  Indeed,  for  a  scheme  in  class  N,  the  interaction  of  weak  waves 
m2 

augments  the  total  variation  norm  at  most  quadratically ,  i.e. 


TVMX  <  TVX  +  ce2(X)  , 

if  as  osc  X  <<  1,  where  the  oscillation  of  X  is  defined  as  the  supremum  of 

K  k 

the  absolute  values  of  wave  magnitudes  6^  and  y  in  X,  and  we  obtain  the 
following  lemma. 

Lemma.  Given  a  scheme  in  class  N  and  a  constant  m  >  0,  there  exists  a 
constant  c(m)  such  that  the  functinal  F1  =  TV  +  c(m)P  satisfies 

F1 (MX)  <  F1 (X) 

if  X  e  rC„  and  if  TVX  «  1. 
m2 

Since  the  functinals  F1  and  TV  are  equivalent  on  states  with  small 
total  variation,  i.e. 

cTVX  <  F^ (X )  <  cTVX 

if  TVX  <<  1,  it  follows  as  a  corollary  that 

k  p 

TVM  X  <  cTVM*X  for  p  <  k  <  q 

if  TVmPx  <<  1  and  if  MkX  lie  in  TC  for  p  <  k  <  q.  Granting  the  lemma, 

m2 

the  problem  of  establishing  stability  in  the  total  variation  becomes  one 

estimating  the  total  variation  norm  along  those  segments  of  the  orbit  which  lie 

„c 


near  simple  waves,  i.e.  in  r  . 

m2 

The  Structure  of  P.  The  potential  P  is  sum  of  n  functionals  Pj,  each 
measuring  the  potential  for  interaction  in  waves  of  a  given  characteristic 
field.  The  functional  Pj  contains  a  constant  coefficient  quadratic  potential 
for  interaction  as  in  the  random  choice  potential  [12]  plus  a  weighted  sum 
corresponding  to  numerical  self-interactions: 

2 


(4.3) 


P  (X)  -  l  aB  +  l  $ ,  ( a,u  )a  . 

j  o  j  a 


Here  ■  Dj(X)  denotes  the  set  of  all  pairs  (a,0)  of  distinct  j-waves  in 
X.  The  weight  A  (a,u  )  depends  on  the  local  base  state  u  through  which  the 
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wave  a  passes,  on  orientation  of  the  normal  of  the  line  segment  which  a 
crosses  and  on  the  way  in  which  the  scheme  under  consideration  approximates  a 
j-simple  wave. 

We  note  that,  in  contrast,  the  random  choice  potential  contains  terms  of 
the  form 

I  Iasi 

A. 

3 

where  A^  denotes  the  set  of  all  pairs  (a, (5)  of  approaching  j-waves  in  a 
given  wave  configuration;  in  the  terminology  of  [12]  a  pair  of  j-waves  is  called 
approaching  if  at  least  member  is  a  shock.  The  presence  of  terms  of  the  form 
(4.4)  is  motivated  by  the  fact  that,  in  the  random  choice  method,  approaching 
j-waves  a  and  6  will  collide  in  a  finite  time  in  the  absence  of  interference 
from  other  waves,  just  as  in  an  exact  solution;  if  a  and  0  are  both  shocks 
then  the  total  variation  typically  increases  at  the  point  of  interaction  by  a 
quantity  on  the  order  |a0|;  on  the  other  hand  if  only  one  wave  is  a  shock,  the 
total  variation  is  reduced  by 

-C(a,0)  +  O(a0)  , 

where  the  cancellation  between  a  and  0  is  defined  by 

C(a,B)  =  2  min(|a|,|0|)  if  sgn  a  *  sgn  0 
=0,  otherwise  . 

cf  [12,  13].  In  constrast,  pairs  of  j-rarefaction  waves  in  conservative 
difference  approximations  can  also  interact  through  numerical  errors  and  one  may 
expect  their  products  to  appear  in  potentials  for  wave  interaction.  Lastly,  we 
remark  that  for  class  N  schemes,  one  has  the  option  of  working  with  potential 
functionals  in  which  the  leading  term  of  (4.3)  is  replaced  by 

I  Ia0|  , 

D, 
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but  in  this  case  additional  terms  involving  cancellation  effects  C(a,B )  appear 
on  the  riy.it  hand  side  of  (4.2).  Such  potentials  are  useful  in  the  context  of 
hybridized  scheme;  for  current  purpose  it  suffices  to  discuss  potentials  of  the 
form  (4.3) 

We  shall  now  turn  to  the  problem  of  estimating  the  total  variation  norm 
along  those  segments  of  the  orbit  which  lie  near  simple  waves.  To  this  end  one 
is  led  to  study  the  local  recession  of  waves  after  interaction;  the  motivation 
is  the  following.  We  recall  that  in  the  exact  solution  operator  three  main 
mechanisms  of  stability  are  present  in  the  form  of  the  cancellation  process 
between  shocks  and  rarefaction  waves  of  the  same  field,  in  the  spreading  of 
individual  rarefaction  waves  and  in  the  recession  of  waves  of  different  fields 


after  interaction.  Now,  if  a  local  state  X^  =  (6  ,Y1#s)  lies  close  to  a 
j-simple  wave,  for  example,  in  the  sense  that 


d  (X, )  <  me  IX,) 
j  P  1 


with  p  =  2  and  m  small  or  with  p  >  2  and  osc  X..  small  then  clearly 
sgn  =  sgn  Yj  and  cancellation  is  absent.  Secondly,  one  expects  that  the 

spreading  of  rarefaction  waves  will  only  be  detected  in  the  framework  of  a 


n-parameter  interpolation  after  several  time  steps.  Therefore,  in  studying  the 
behavior  of  a  class  N  scheme  in  two  consecutive  time  steps  acting  on 
configurations  near  simple  waves,  it  is  natural  to  investigate  the  wave 


recession  process.  For  this  purpose  we  shall  construct  a  functional  T  which 


measures  the  potential  for  transverse  wave  interactions  and  satisfies 
(4.5)  T(MX )  -  T (X )  <  -CT(X)  +  cd2(X)  +  ce3(X)  , 

if  TVX  <<  1.  Here  t  records  the  sum  of  products  of  all  incoming  transverse 
waves  in  X : 

T (X )  =  l  {|v*| |v*|  :  j  <  i  and  <  k  <  «}  . 
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In  particular,  if  X  lies  in  T  (for  appropriate  m)  then  T  is  reduced  in 

m3 

one  time  step  by  a  quantity  on  the  order  of  x  modulo  the  approximate  rate  of 
entropy  production,  i.e. 

T(MX)  -  T(X)  <  -cx(X)  +  ce3(X)  , 

if  TVX  <<  1  and  X  e  Thus  if  X  e  I*  the  functional  T  experiences 

virtually  the  same  reduction,  modulo  e^,  as  the  exact  solution  operator.  In 

this  regard  we  note  that  x  includes  products  of  both  approaching  and  receding 

waves  where,  in  the  standard  terminology,  a  pair  (6.,y„)  is  called  approaching 

j  * 

if  6.  lies  to  the  left  of  y  and  either  j  >  i  or  j  =  l  and  at  least  one 
J  * 

member  is  a  shock  and  receding  if  6 .  lies  to  the  left  of  y.  and  j  <  l.  Now 

j  x 

if  a  local  state  is  "close"  to  the  simple  waves  then  x  is  in  fact  on  the  order 
of  the  sum  of  just  the  products  of  those  incoming  waves  which  are  approaching. 
The  latter  quantity  is  precisely  the  amount  by  which  the  potential  for  the 
random  choice  method  is  reduced  in  a  local  interaction  cf.  [12 J.  We  conclude 
that  for  schemes  in  class  N,  a  potential  of  the  form  Q  =  cP  +  cT  satisfies 
(4.6)  Q(MX)  -  Q(X)  <  -cx(X)  -  cd2(X)  +  ce3(X) 

if  TVX  <<  1. 

The  form  of  the  right  hand  side  of  (4.6)  leads  one  to  study  the  action  of 

class  N  schemes  on  the  following  sets  of  states  W_  contained  in  T  , 

»P  mp 

W  *  {X  :  x(X)  +  d2(X)  <  me  (X)}  . 
mp  p 

If  p  >  2  or  if  p  =  2  and  m  is  small,  one  can  regard  a  point  X  in  wmp 

as  a  weakly  interacting  state  in  the  sense  that  the  total  amount  of  transverse 

wave  interaction  in  X  plus  its  distance  squared  to  the  corresponding  simple 

waves  is  relatively  small.  In  the  case  p  =  2  we  obtain  a  bound  on  the  total 

lc  c 

variation  along  those  segments  of  the  orbit  MX  which  lie  in  W 


Lemma.  Given  a  scheme  in  class  N  and  a  constant  m  >  0  there  exist 
constants  c^  and  Cj  depending  on  m  such  that  the  functional 
F  *  TV  +  c^P  +  c^T  satisfies 

(4.7)  F{ MX )  <  F(X ) 

if  TVX  <<  1  and  if  X  e  w°  . 

m2 

Again,  it  follows  as  an  immediate  corollary  that 

(4.8)  TVMkX  <  cTVMPX  for  p  <  k  <  q 

if  tvmPx  <<  1  and  if  MkX  lies  in  wc„  for  p  <  k  <  q. 

m2 

In  Section  6  we  shall  describe  a  procedure  for  perturbing  the  generating 
function  of  class  N  schemes  in  such  a  way  as  to  damp  the  numerical  reflected 
waves  which  are  produced  by  the  interactions  of  configurations  near  simple 
waves.  For  the  perturbed  schemes  estimates  of  the  form  (4.7)  and  (4.8)  hold 

Jr 

provided  the  X  and  M  X  respectively  lie  in  the  complement  of  the  much 
smaller  third  order  states  Wm3 .  As  we  remarked  above  the  sets  Wm3  are  of 
particular  interest  since  e3(X)  represents  the  approximate  rate  of  entropy 
production. 

Structure  of  T.  The  functional  T  consists  of  two  terms.  The  first  represents 
the  standard  quadratic  potential  for  approaching  waves  of  different 
characteristic  fields  and  the  second  represents  a  numerical  potential  for  self¬ 
interactions  within  groups  of  waves  associated  with  pairs  of  adjacent  mesh 
points: 

T(X )  *  l  |a6l  +  l  ^(a,B,u  ,ufl)|a0|  . 

A  0  6 

Here  A  =  A(X)  denotes  the  set  of  all  pairs  (a,0)  of  approaching  waves  of 
different  fields  and  the  weight  i|>  has  the  following  property:  if  a  is  a 
j-wave  and  0  a  k-wave,  k  *  j,  crossing  the  west-south  side  or  the  south-east 
side  of  a  mesh  diamond  with  base  state  ug  respectively  then 


*  =  pjk<V  °r  *  -  VV 

for  certain  smooth  functions  Pjk  and  q^,  which  characterize  the  numerical 
self-interactions  among  the  waves  of  the  associated  Riemann  problems,  otherwise 
i|i  =  0.  The  details  of  the  construction  of  T  are  given  in  the  next  section. 

5.  POTENTIAL  FUNCTIONALS. 

In  this  section  we  shall  construct  the  functionals  Pj  and  T  discussed 

in  Section  4.  The  construction  of  Pj  is  based  on  the  following  fact. 

2 

Consider  the  linear  mapping  of  R  defined  by 

(5.1)  a  =  Av?  A  =  1  -  w  t 

p  1  -  T 

where  v  =  (6,y)  and  a  =  (a,0).  Such  mappings  arise  from  the  leading  term  of 

(3.1)  with  a  fixed  value  of  the  local  base  state  ug.  A  simple  computation 
shows  that  quadratic  form 

Q ( v)  =  ay2  +  b52  +  y5 

satisfies 

(5.2)  Q(SAV)  -  Q(AV)  =  -(y  +  t)(1  -  y)(b  -  |)(5  -  0y)2 
if  the  constants  a  and  b  satisfy 

(5.3)  a  -  1/2  =  0(b  -  1/2);  0  =  (1  -  t)/(1  -  u)  . 

Here, 


Indeed,  if  the  left  hand  side  of  (5.2)  maintains  one  sign  for  all  v  then  (5.2) 

and  (5.3)  necessarily  hold.  In  the  context  of  schemes  in  class  N  this 

observation  yields  the  following  result.  Fix  the  local  base  sate  ug  and  put 

A  =  A , (u  ),  u  =  U,(u  )  and  T  =  t  (u  )  . 

3  s  j  s  j  s 

It  follows  from  (3.7)  and  (3.8)  that 
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(1  -  t.)/(1  -  U.)  =  9 , (u  )  i  {a  +  a  A  (u  )}/{B  +  B  A.(u  )} 

]  3  3  s  txjs  t  x  3  s 


Thus,  if  we  set 


Ma,u  )  *  b ,  (u  ) 
3  a  3  s 


if  a  is  a  j-wave  crossing  the  ws  side  of  a  mesh  diamond  with  base  state 
ug  and 

$  .  (a,u  )  =  a,(u  ) 

3  a  3  s 

if  a  is  a  j-wave  crossing  the  se  side  of  a  mesh  diamond  with  base  state 

n 

ug  then  the  functional  P  =  £  P  with  P.  defined  by  (4.3)  satisfies  the 

j  =  1  3 

inquality  (4.2)  if  and  b^  satisfy 

(5.4) 

and  if  TVX  <<  1.  we  note  that  changes  in  the  arguments  of  the  weights 


(a..  -  1/2)  =  0..(b.  -  1/2);  b.  >  1/2  , 


leads  only  to  cubic  contributions  which  can  be  absorbed  into  e3(X).  We  note 

that  the  restriction  (5.4)  depends  only  on  the  eigenvalues  A  and  the  choice 

of  normals  a  and  3  and  leads  to  a  one-parameter  family  of  potentials  P^. 

This  is  in  accordance  with  the  fact  that  any  finite  difference  scheme  with 

smooth  generating  function  consistent  with  the  equations  is  at  least  first  order 

accurate  on  smooth  solutions:  if  a  mesh  diamond  is  embedded  in  an  exact  simple 

wave  and  if  the  values  of  the  exact  solution  at  the  west,  south  and  east 

vertices  are  substituted  into  the  generating  function  <p,  then  the  value 

produced  by  $  differs  from  the  value  of  the  exact  solution  at  the  north  vertex 

2 

by  a  quantity  on  the  order  of  Ax  . 

The  potential  T(X)  for  transverse  interactions  is  constructed  as 

follows.  We  begin  by  considering  a  slightly  more  general  potential  of  the  form 

(5.3)  T(X )  =  s  l  |aBl  +  w  l  ! aB I  +  \  p(a,B)|aB|  +  £  r(a,B)|aBl  . 

A  R 


Here  A(X)  and  R(X)  denote  respectively  the  sets  of  all  approaching  and 
receding  pairs  of  waves  in  X  such  that  each  member  wave  is  associated  with  a 
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different  pair  of  adjacent  mesh  points;  p(ct,B)  =  p^k(ug)  if  a  is  a  j-wave 
and  8  a  k-wave  crossing  the  ws  side  of  a  mesh  diamond  with  base  state  ug, 
p  =  0  otherwise;  r(a,8)  =  rj^us)  if  a  is  a  j-wave  and  8  a  k-wave 
crossing  the  se  side  of  a  mesh  diamond  with  base  state  ug,  r  =  0  otherwise. 
We  shall  show  that  T(X)  satisfies  (4.5)  for  an  appropriate  choice  of 
coefficients,  in  particular,  for  choices  such  that  s  is  constant,  w  vanishes 
and  p^  and  rjk  are  smooth  positive  functions  of  ug.  To  this  end  we  write 
the  linear  map  (5.1)  in  the  form 

a  =  y  +  a(6  -  0y) 

8=6-  a ( 6  -  0y),  a  =  1  -  y  , 

in  order  to  display  the  deviation  for  the  corresponding  simple  wave,  6  =  0y. 
Consider  a  discrete  wave  interaction  associated  with  a  mesh  diamond  having  a 
local  base  state  ug.  Fixing  the  value  ug,  the  incoming  and  outgoing  waves 
v  =  (6,y)  and  a  =*  (a, 8)  are  related  as  follows  modulo  a  term  on  the  order 


|v| 


2 


(5.4) 


+  V  *1  5  Vsi  -  W  - 

6.  -  -  *. 


where  a,  =  a.(u„)  and  0J=0,(u).  In  order  to  establish  (4.5)  it  is 
J  J  s  j  j  s 

sufficient  to  show  that  the  following  incoming  and  outgoing  potentials 
associated  with  a  single  discrete  wave  interaction. 


Vv> '  ■  S  Vk  *  * 1  s)Tk  * 1  p)kVk  *  J  r)kVk  ■ 
V0>  -  •  l  *  “  l  *  I  ♦  I  r-nVk  ' 


with  summations  taken  over  indices  j  <  k,  satisfy 


(5.7) 


T  (0)  -  Ti(v) 


T(V)  d=f  I 
j<k 


<  -const  t(v)  +  const  £  |t|/  i|/  | 

j<k  3  * 

,VjMV'  ,Vt’  "  ’V  +  ’V 


provided  that  (5.4)  holds  and  the  coefficients  are  properly  selected.  We  note 
that  the  quadratic  correction  terms  to  (5.4)  as  well  the  changes  in  the  local 
base  state  introduce  only  admissible  cubic  terms  into  formulas  (5.7).  Lastly, 
we  remark  that  without  loss  of  generality  one  can  restrict  attention  to  the  case 
where  6^  and  y^  are  non-negative  for  j  =  1,2,...,n  since  the  general  case 
follows  by  replacing  <5^  and  y  ^  with  |6^|  and  |y^l  respectively. 

Next,  we  shall  describe  the  calculations  which  lead  to  the  restrictions 
on  a'w'Pjk  and  rjfc  which  guarantee  (5.5).  In  the  following  all  summations 
are  taken  over  indices  j  <  k.  Substituion  of  (5.4)  into  (5.6)  yields 
T0(a)  -  T± (v)  -  (s  -  w)  l  6jYk  -  Y^  +  l  (rjk  +  Pjk  -  s  -  w)^^ 

+  l  (S°j  -  Pjk6j  -  +  rjkV\  +  Z  (w6k  -  PjA  *  ^k  +  rjkV*j 

A  brief  calculation  show  that 


"  *  V"  *  V'V*  ■  Vk>  ’  (9j  -  V'S  +  V‘5k  +  V 


*  11  +  V<\  *  Wa)  ■  "  *  V,T)  *  VVs 


and  therefore 


T. (o)  -  T. (v)  **  (s-w)  l  (0.  -  0  )(6.  +  Y.  )(«. 
0  i  ifcjjk 

(5.8)  +  l  [{w  -  p  +  (s-w)/a  (1+0  )}<$  +  {r 

jK  J  J  *  J* 

+  I  {s  -  pjk  +  (s-w)/ak(1+0k)}6j  +  (rjk 


The  condition  that  the  coefficients  of  and 

summations  of  (5.8)  vanish  when  6^  =  0  y^  and 
to  the  following  pair  of  linear  equations  for  r 
(specified)  values  of  s  and  w 


+  Y.  )/( 1+9 , ) ( 1+6  ) 
k  j  k 

-  s  +  (s-w) /a  ( 1+0  )  }y.1  <P  . 

j  j  K  1 

-  w  -  ( s-w) /a^ ( 1+9k ) >Y ^ 

in  the  second  and  third 
6k  =  ®kYk  respectively  leads 
jk  and  p^k  in  terms  of 
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;T 


where  r =  (s  -  w)/a.(1  +  9.).  The  system  (5.9)  has  a  unique  solution  since 
Ik  3  D 

0  .  *  0  if  j  *  k.  In  particular  we  note  that  if  we  take  w  =  0  and  s  to 
3  k 

be  a  positive  number  then  the  solutions  r^  and  p^  are  positive  and  we 
obtain  an  identity  of  the  form 

T  (a )  -  T.  (v)  =  si  (0.  -  0.  )(6.  +  y.)(6,  +  Y.  )/(1  +  9.)(1  +  0,  ) 

0  i  jkjjkk  j  k 

+  l  [s{l/a.(1  +  0.)  -  1 /a  (1  +  6  )  +  1}  -  2p  ]*.* , 

1  1  *  K  1  k 

in  which  the  first  summation  involves  a  negative  coefficient  and  the  second  is 
on  the  order  of  the  square  of  the  distance  to  the  corresponding  simple  waves. 

It  follows  that  the  functional  T(X)  of  the  form  (5.3)  satisfies  the  desired 
estimate  (4.5)  if  we  take  w  =  0,  s  positive  and  r^k  and  p^  as  the 
solutions  of  (5.9).  Here  the  functions  r ^  and  p^k  will  depend  smoothly  on 
the  local  base  state  ug  since  the  equations  (5.9)  depend  smoothly  on  ug.  We 
note  that  an  analogous  functions  with  w  *  0  can  be  constructed  but,  as  we 
shall  see,  such  a  functional  is  less  convenient  for  the  purposes  of 
hybridization  with  the  random  choice  method.  Indeed,  the  act  of  attributing  a 
potential  for  interacting  of  the  form  w|aP|  to  a  pair  of  receding  waves 
a  and  0  in  the  random  choice  method  leads  to  several  unnecessary  terms  which 
are  awkward  to  handle.  These  terms  are  simply  avoided  by  setting  w  =  0. 

6.  REFLECTION  AND  CANCELLATION  OF  WAVES. 

We  observed  in  Section  4  that  a  functional  of  the  form  TV  +  c^  +  c^T  is 
decreasing  along  these  segments  of  the  orbit  M  X  with  small  total  variation 
which  lie  in  the  complement  of  the  weakly  interacting  states  wm2*  In  this 


section  we  shall  construct,  for  a  given  generator  $  in  class  N,  a  quadratic 
perturbation 

2  2 

q  =  q(w,s,e)  =  0{ | w  -  s|  +  Is  -  e|  }  , 

with  the  property  that  the  reflected  waves  produced  by  the  generator  $  +  q  on 
local  states  near  simple  waves  are  third  order  with  respect  to  the  incoming 
magnitudes  and  the  property  that  the  marching  map  M  associated  with  <j>  +  q 
satisfies 

2 

TVMX  -  TVX  <  const. {t(X)  +  d  (X)  +  e3(X)}  , 
if  osc  X  is  small,  i.e.  in  one  time-step  the  total  variation  norm  can  not 
increase  by  more  than  a  quantity  on  the  order  of  the  corresponding  reduction  in 
the  potential  Q  =  c.jP  +  C2T'  modulo  a  quantity  on  the  order  e3(X)  of  the 
approximate  rate  of  entropy  production.  It  follows  that  for  appropriate  m  the 
functional  F  =TV  +  Q  is  decreasing  along  those  segments  of  the  orbit  with 
small  total  variation  which  lie  in  the  complement  of  Wm3  and  that 

F( MX )  -  F(X )  <  const. e3(X) 

for  arbitrary  configurations  X  with  small  total  variation. 

In  this  section  we  shall  also  describe  the  effective  cancellation  between 
shocks  and  rarefaction  waves  of  the  same  field  which  exists  for  class  N 
schemes  and  their  perturbations  and  compare  it  with  the  corresponding 
cancellation  occurring  in  the  exact  solution.  To  this  end  we  shall  begir.  by 
recalling  that  for  schemes  in  N  the  outgoing  and  incoming  magnitudes  of  a 
local  interaction  and  related  by  formulas  of  the  form 


(6.1) 


a.  =  (1  -  U  . ) 6  .  +  t.y.  +  p.(v);  p.(v) 

1  11113  3 

8.  =  y.5.  +  (1  -  T.)Y.  +  q  .  (  V )  ;  q  .  ( V ) 
113  3  11  1 


0( | V | 2 )  , 

0(|v|2)  , 


where  v  =  (v  ,...,v  )  and  the  functions  y.,T.,p.  and  q  depend  smoothly  on 
In  113  1 

the  base  state  of  the  associated  mesh  diamond.  The  process  connecting  the 
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incoming  and  outgoing  magnitudes  =  (6^,y,)  and  a =  (a.,S.)  can  be  viewed 


as  a  composition  of  binary  interactions  by  writting  (6.1)  in  the  form 


a  =  a  (v)  =  l  a . ( v  )  +  0{t(v)} 
3  3  k=1  3  k 


(6.2) 


n 


8.  =  8.(v)  =  T  B.(v  )  +  o{t(v)}  . 
3  3  k=1  3  k 


Here  e  Rn  represents  the  vector  whose  kth  component  equals  and  Whose 


jth  component  j  *  k  vanishes;  the  local  base  state  is  regarded  as  fixed  and 
t ( v )  measures  the  total  transverse  wave  interaction: 


T(V)  *  l  I V .  |  | V  |  . 

j<k  3  * 


Indeed,  it  is  easy  to  show  that  a  formula  of  the  form  (6.2)  holds  for  an 
arbitrary  function  a  **  a(v)  vanishing  at  the  origin  and  having  a  bounded 
second  derivative. 

We  shall  construct  the  perturbation  q  in  such  a  way  that  incoming 
k-waves,  produce  outgoing  j-waves,  j  *  k,  from  the  generator  <j>  +  q  which  are 
on  order  of  the  distance  squared  to  the  corresponding  k-simple  waves  modulo  a 
cubic  error  and  in  such  a  way  that  conservation  of  wave  magnitudes  holds  in  the 
jth  field  modulo  a  cubic  error,  i.e. 

Iaj(vk)|  +  I8j(vk)l  <  const.  |  6  ^  -  ek^us)3'k|2  +  const.  I  |  3  , 

A  A  n 

aj<Vj)  +  03(Vj)  =  5j  +  Yj  +  0( I v j |  )  . 

It  follows  that  (6.2)  can  be  written  in  the  convenient  form 

A 

=  (1  -  ^><5.  +  TjYj  +  Pj<vj)  +  °{g<X)} 

A 

0,  =  M  j5  ,  +  (1  -  T  .  )  Y  .  +  qjv.)  +  0{g(X ) } 


A  A  A  . 

p.(v.)  +  qj(Vj)  =  0< 'Vjl  ’ 

g(X)  =  T ( X )  +  d2(X)  +  | v | 3 

and  X  =  (<S,Y,u  )  e  R3"  specifies  the  incoming  configuration.  One  can  then 
s 

employ  (6.3)  to  display  the  effective  cancellation  experienced  by  interacting 
j-shocks  and  j-rare faction  waves  by  noting  that 

la.  I  +  16.1  <  |6.l  +  (Yj  I  -  C^v)  +  0{g(X ) } 

where 

C.(v)  =  C( 6  ,y . )  -  C(a.,B.) 

3  j  3  3  3 

{2min( |x| , |y | )  if  sgn  x  *  sgn  y 
0  otherwise 

Furthermore,  it  is  not  difficult  to  show  that  the  effective  cancellation  C.(v) 

3 

in  the  jth  field  is  bounded  below  by  a  fraction  of  the  cancellation  C(6^,y.) 

which  occurs  in  the  interaction  of  waves  6_.  and  Y ^  in  the  random  choice 

method,  modulo  a  quadratic  term: 

C.(v)  >  k.C(<S.,YJ  -  0{g(X ) } 

3  3  3  j 

kj  =  ^  min(  W  ^  ,  1  -  V  j » t  ^  ,  1  -  t_.)  . 

Thus,  we  obtain  the  following  estimates  on  the  outgoing  waves  of  local  and 
global  interactions: 

|a  I  +  |0  |  <  |6  |  +  | Y . !  ~  k.C(6  ,Y.)  +  const. g{(S,Y,u  )} 
3333333  s 

TVMX  -  TVX  <  -const. C(X)  +  const. (d2(X)  +  T(X)  +  e3(X)} 

k  k 

where  the  cancellation  in  a  global  configuration  X  =  {(5  ,y  ,sk)}  is  defined 
by 

C(x)  =  l  C(6k,Y*)  . 
j,k  3  1 
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In  contrast,  the  interaction  of  a  j-shock  and  j-raref action  wave  in  the  random 
choice  method  leads  to  the  absorption  of  the  smaller  wave  by  the  larger  and  a 
concurrent  decrease  in  the  total  variation  norm  equal  to  twice  the  magnitude  of 
the  smaller  wave  modulo  quadratic  terms:  the  outgoing  j-wave  is  related  to 

the  incoming  j-waves  S_.  and  Y_.  by  a  formula  of  the  form 

I  e  J  =  !<S^I  +  I Y  j  I  -  C(<S_.,y..)  +  0(  |  v  (  2  ) 
cf.  [13].  For  conservative  difference  schemes  the  presence  of  fractional 
cancellation  leads  to  the  persistance  of  small  oscillations,  corresponding  to 
alternating  sequences  of  shocks  and  rarefactions,  over  several  time  steps. 

We  shall  next  construct  an  appropriate  perturbating  q  for  a  given 
generator  <{>  in  class  N.  Consider  a  discrete  interaction  associated  with  a 


mesh  diamond  and  regard  the  base  state  ug  as  fixed.  Suppose  that  the  incoming 
waves  consist  of  just  a  pair  of  j-raref action  waves,  i.e.  suppose 


5.  >  0,  y ,  >0  and  V  =  0,  k  *  j.  Here  the  values  at  the  west  and  south 
ID  * 

vertices  u„  and  u_  lie  on  the  j-wave  curve  r.(s)  through  u  ,  the 

W  o  j  S 

integral  curve  in  Rn  of  the  right  eigenvector  r^(u),  and  satisfy 

X.(u  )  <  X  fu  )  <  X  (u  )  . 

3  w  j  s  j  e 

If  the  generator  <f>  were  to  produce  a  point  on  then  reflected  waves  are 

absent: 


Oj  +  3.  =  5.  +  Yj;  ak,0k  =  O'  k  *  j  * 

However,  in  general,  the  point  n_,  =  n_.(5^,Y^,s)  produced  by  (j>  for 

configurations  with  =  0,  k  *  j  lies  within  a  distance  on  the  order  of 
2 

|v.|  from  the  nearest  point  of  r_j(s)  which  we  shall  denote  by 

We  shall  first  construct  a  perturbation  q.  which  reduces  to  m.  -  n.  if  v 

3  J  J  j 


forms  a  discrete  j-rarefaction  wave  and  vanishes  if 


|5  -  8  (u  )Y , I  >  const. ( 1 6 . |  +  |  y  - 1  >  • 

j  3  s  j  3  3 

Fix  m  >  0,  let  z(x)  be  a  smooth  even  function  equal  to  one  of  x  <  m  and 

equal  to  zero  if  x  >  2m  and  put 

q.(v)  =  q.(v  )  =  Z{(5.  -  6  .  (u  )Y. )/(<$.  +  Y  , ) )  (m .  -  n.}  . 

3  3  3  3  3  s  3  3  3  3  3 

The  desired  perturbation  q  is  obtained  by  taking  m  small  and  defining 

(6.4)  q  =  I  d . ( v)q . ( v) j  d.  =  II  z{(62  +  y*)/(62  +  y2))  • 

3=1  3  3  3  *  *  3  3 

We  note  that  while  the  formula  (6.4)  for  q  is  motivated  by  the  case  where  the 
incoming  configuration  consists  of  3ust  a  pair  of  incoming  j-raref action  waves 
it  has  the  desired  effect  for  general  configurations.  In  particular  we  note 
that  the  support  of  q  is  contained  within  the  set  of  all  states 
(6,Y»ug)  e  R3n  with  the  property  that  there  exists  an  index,  say  j,  such 
that 

|vk!  4  4m|vj|,  k  *  3 

|6  -  8  (u  )Y . I  <  2m( 1 6 . |  +  | Y  .  I ) ?  sgn  6.  =  sgn  y.  • 

3  3  s  3  33  3  3 

Furthermore  a  straight  forward  calculation  shows  that  q  vanishes  together  with 

its  first  derivative  in  v  at  v  =  0  and  has  a  bounded  second  derivative: 

2 

roughly  speaking  q  makes  a  change  on  the  order  of  e  over  a  distance  on  the 
order  of  e . 

7.  THE  RANDOM  CHOICE  SCHEME. 

We  shall  briefly  describe  the  generating  function  of  the  random  choice 
method  together  and  the  potential  functionals  which  are  used  to  obtain  a  uniform 
estimate  on  the  total  variation  norm  of  the  corresponding  difference 
approximations.  We  shall  compare  these  functionals  with  the  functionals 
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constructed  in  Section  5  for  schemes  in  class  N  and  describe  certain 
modifications  of  the  random  choice  generating  function  which  facilitate  a 
hybridization  with  schemes  in  class  N. 

The  random  choice  method  can  be  based  on  a  centered  diamond-shaped  stencil 

as  follows.  Let  R(x/t;u-,u+)  denote  the  classical  solution  of  the  Riemann 

problem  with  data  (u",u+)  [20],  and  let  (y^)  be  a  sequence  equidistributed 

in  the  interval  (-1,1).  Consider  a  grid  whose  stencil  has  vertices  of  the  form 

n  =  {x,(m  +  1)At}  w  =  (x  -  Ax,mAt} 

s  =  {x,(m  -  1)At}  e  =  {x  +  Ax,mAt} 

where  x  is  an  integral  multiple  of  Ax.  The  value  of  the  random  choice 

approximation  at  the  north  vertex  depends  only  on  the  values  at  the  west  and 

east  vertices  and  the  corresponding  element  ym;  it  is  obtained  by  solving  the 

Riemann  problem  with  data  (uw,ufi)  and  sampling  the  value  of  the  solution  at 

time  t  *  At  and  position  x  *  y  Ax,  i.e. 

m 

U  =  u(u,u,y)  =  R(y  Ax/At ;u,u)  , 
n  nwem  m  we 

the  C-F-L  condition  is  enforced  in  the  standard  fashion.  For  convenience  one 
can  associate  with  each  grid  function  produced  by  the  randan  choice  method  a 
piecewise  constant  function  which  assumes,  on  each  small  parallelogram  centered 
at  a  mesh  point,  the  corresponding  value  of  the  grid  function.  We  note  that  the 
random  choice  method  in  its  original  formulation  [12]  involves  approximate 
solutions  which  are  exact  in  strips  of  the  form  mAt  <  t  <  (m  +  1)At.  However, 
it  is  not  difficult  to  show  that  if  a  sequence  of  such  approximate  solutions  is 
stable  in  the  total  variation  norm  and  convergent  pointwise  a.e.  then  the 
corresponding  piecewise  constant  functions  described  above  converge  pointwise 
a.e.  to  the  same  limit. 

In  order  to  discuss  the  functionals  introduced  by  Glimm  for  the  random 
choice  method,  we  shall  briefly  recall  the  laws  governing  the  local  wave 


interactions  therein.  We  first  note  that  the  outgoing  wave  magnitudes  of  a 

local  interaction  associated  with  a  single  mesh  diamond  have  the  following 

structure.  Let  e  =  e.(u  ,u  )  denote  the  magnitudes  of  the  waves  in  the 
j  j  w  e 

Riemann  problem  with  data  (uw,ue).  There  exists  an  index  m  such  that 


(7.2) 


=  e j ,  3^  *>  0  ;  j  ■  1,...,m  -  1 

aj  =  °  =  Gj?  J  =  m  +  1'***'n 


Furthermore,  if  e  <  0  then  either  a  =  e  and  0  =  0  or  a  =  0  and 

m  m  m  m  m 

8  *  e  accordingly  as  the  sample  point  y  Ax/ At  lies  to  the  right  or  left  of 

ram  m 

the  m-shock  of  the  Riemann  problem.  If  e  >  0  then  either  the  aforementioned 

m 

relation  holds  or  the  sample  point  splits  the  m-rarefaction  wave  of  the  Riemann 

problem  and  produces  waves  a  >  0,  8  >0  which  satisfy  a  +  0  »  e  .  In 

m  m  m  m  m 

contrast  to  conservative  difference  schemes,  local  interactions  in  the  random 
choice  method  can  only  increase  the  number  of  waves  (by  one)  if  a  rarefaction 
wave  is  split  and  such  splitting  is  not  accompanies  by  wave  amplification  in  the 
sense  that 


U .  I  =  15.1  +  lY.I  -  C(<5  ,y  )  +  0{d(5,y)} 

(7.3)  3  3  3  33 

D(6 ,y)  -  l  :  \  and  Yj  approach}  , 

in  particular,  two  rarefaction  waves  of  the  same  field  do  not  approach,  cf. 

Section  5.  It  then  follows  from  (7.3),  together  with  the  local  recession  of 

waves  after  interaction  as  expressed  by  (7.2),  that  the  functional 

Q(X)  =  J  { | ctB  1  :  a, 8  approaching  in  x} 

compensates  for  wave  amplification,  i.e. 

Q( MX )  -  Q(X )  <  -const.  I  D(6k,6k) 

def 

F( MX )  -  F(X )  <  0;  F  =  TV  +  const.  Q  , 


k  k 

if  X  =  {(6  ,y  ,s^)}  has  smaH  total  variation  [12].  We  note  that  F(x)  is 
independent  of  the  base  states  of  X,  their  influence  is  absorbed  into  the  last 
term  on  the  right  hand  side  of  (7.3). 

The  local  recession  of  waves  also  plays  a  central  role  in  our  hybridized 
methods.  In  this  connection  we  shall  next  remark  on  several  ways  in  which  the 
structure  of  the  equidistributed  sequence  (y^)  influences  the  geometry  of 
interacting  waves  and  we  shall  present  two  modif ications  of  the  random  choice 
generator  which  guarantees  that  the  local  recession  of  waves  occurs  in  a  minimal 
number  of  time  steps.  For  simplicity  let  us  restrict  our  attention  to  systems 
of  two  equations.  Consider  a  corresponding  exact  solution  which  consists  of  two 
interacting  shocks  of  different  fields,  cf.  the  solution  u,  discussed  in  the 
subsection  on  self-interactions  in  Section  3.  For  a  solution  of  the  form  u, 
the  random  choice  approximations  consist  of  exactly  two  shocks  in  each  strip 
mAt  <  t  <  (m  +  1)At?  for  a  typical  equidistributed  sequence  the  shocks 
approach  during  a  time  interval  of  the  form  (0,pAt),  interact  at  t  =  pAt, 
remain  adjacent  in  (pAt,qAt)  and  separate  in  (qAt,°°).  The  shocks  remain 
adjacent  in  (pAt,qAt)  if  the  sample  point  lies  strictly  to  the  right  or  left 
of  the  waves  in  the  associated  Riemann  problem;  of  course,  since  the  sequence  is 
equidistributed  this  can  not  happen  for  an  arbitrarily  long  string  of  its 
elements,  hence  the  waves  eventually  separate. 

Certain  technical  problems  associated  with  delayed  wave  recession  in 
hybridized  schemes  can  be  avoided  by  using  an  equidistributed  sequence  such  that 
at  least  one  element  in  every  string  of  m  consecutive  elements  corresponds  to 
a  point  in  the  wake  region  between  characteristic  fields.  In  order  to  describe 
the  construction  of  such  sequences  let  us  restrict  our  attention,  for 
concreteness,  to  systems  with  symmetric  eigenvalues,  i.e.  X^Cu)  =  -X2<u);  the 
general  case  is  handled  analogously.  We  note  that  systems  of  the  form  (1.4) 
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have  symmetric  eigenvalues.  Let  us  suppose  that  all  analysis  takes  place  in  the 
neighborhood  N(u)  of  a  fixed  vector  u  e  Rn  and  write  the  interval 
I  =  (-1,1)  in  the  form 

I=I1UI2UI3 

where  1^  =  ( — 1 , — 1 /2 ] ,  I2  =  (-1/2, 1/2)  and  Ij  =  [1/2,1).  Consider  an 
equidistributed  sequence  (y^)  in  I  which  satisfies 

(7.4)  y2k  e  i2  for  all  k;  y2k+1  e  I.,  for  k  even;  y  2k+-|  e  I  3  for  k  odd. 

Such  sequences  can  be  constructed,  for  example,  by  starting  with  a  sequence 
zk  equidistributed  in  (0,1),  scaling  and  translating  in  the  obvious  way  tc 
obtain  three  sequences  (a^),  ( b^ )  and  (c^)  equidistributed  in  1^,  I2, 
and  I3  and  then  arranging  their  elements  in  the  alternating  manner  indicated 
by  (7.4),  i.e. 

Y4k  '  bk;  y4k+1  “  %'  Y4k+2  =  bk+l'  y4k+3  =  *k  ' 

Now,  if  the  C-F-L  number  is  choser;  so  that 

{X^ (y) At/Ax  ;  u  e  N(u)}  C  I1  and  {X2(u)At/Ax  :  ue  N(u)}  1 ^ 

then  every  second  member  of  the  sequence  (y^)  corresponds  to  a  sample  point 

y^At/Ax  which  lies  in  the  wake  region  in  the  sense  that 

max  X.  <  y,  At/Ax  <  min  X„  . 

N(u)  1  k  N(u)  2 

In  working  in  the  context  of  general  systems  or  in  the  absence  of  any 
restrictions  on  the  C-F-L  number  (except  that  it  be  less  than  one)  one  will 
obtain  an  equidistributed  sequence  such  that  every  mth  element  m  >  2  lies  in 
the  wake  region. 

We  note  that  the  process  above  which  construct  yk  from  z^  does  not 
alter  the  error  of  approximation  of  z^  in  the  following  sense.  If  the  average 
number  of  elements  of  z^  with  index  <  n  in  a  given  subinterval  J 


A 


approximates  half  the  length  of  J  within  an  error  on  the  order  of  ij>(n),  i.e. 

I S{n, j; (z^) }/n  -  |  J|  /2  |  <  c(J)iMn) 

where,  without  loss  of  generality  'J'(n)  vanishes  monotonically ,  then  the  same 
rate  of  approximate  obtains  for  (y^),  i.e. 

|s{n,J;  (yk)}/n  -  |J|/2|  <  d(J)iHn) 

for  an  appropriate  constant  d(J)  depending  only  interval  J.  For  a  typical 

-1  /2 

equidistributed  sequence  one  has  ip(n)  =  n  ,  while  quadratic  irrationalities 
lead  to  sequences  for  which  ip(n)  =  log  n/n.  For  technical  reasons,  concerning 
hybridized  schemes,  it  will  be  necessary  for  us  to  employ  equidistributed 
sequences  with  a  rate  i|i(n)  satisfying 
(7.5)  lim  iji(n)n3/4  =  0 

n+<» 

in  order  to  show  that  limit  of  the  associated  difference  approximations  is  an 

exact  solution.  As  we  shall  see  in  Section  11  the  restrictin  (7.5)  is  not 

necessary  for  the  stability  of  the  hybridized  schemes. 

Next,  we  shall  describe  a  second  modification  of  the  random  choice 

generator  which  has  the  desirable  effect  of  emphasizing  the  local  recession  of 

waves  after  interaction.  For  this  purpose,  let  us  consider  a  random  choice 

interaction  for  a  system  of  two  equations  and  suppose  that  the  sample  point 

y, Ax/At, At)  does  not  lie  between  the  waves  e.  in  the  solution  of  the  Riemann 
k  3 

problem  with  data  <uw»ue)  obtained  from  the  west  and  east  vertices  of  the 
associated  mesh  diamond.  In  this  situation  one  of  four  configurations  obtains 
for  the  outgoing  waves  (ct,8):  both  cross  the  wn  side;  both  cross  the  ne 
side;  is  a  rarefaction  wave  split  by  y^,  i.e., 

al  >  °'  “2  =  °'  B1  >  °'  °1  +  ^  =  G1  ; 
e2  is  a  rarefaction  wave  split  by  y^ ,  i.e.. 


a2  >  °'  62  >  °'  S1  “  °'  °2  +  6  2  =  £2  ’ 
If  any  of  the  above  four  configurations  occur  and  if 


(7.6) 


1  e1  e2 1  >  eo  »  U1  I  +  le2t  <  La 


with  e  and  L  given  positive  constants,  we  replace  the  original  element  y^ 
by  zero.  Thus,  in  the  case  of  outgoing  waves  whose  magnitudes  are,  so  to  speak, 
of  order  a  we  enforce  the  immediate  recession  of  waves  if  it  does  not  occur  in 
the  original  choice  of  equidistributed  sequence. 

The  effect  of  such  a  modification  on  the  error  associated  with  the  random 
choice  difference  approximations  can  be  analyzed  as  follows.  Let  £)(i,j) 
denote  the  mesh  diamond  centered  at  (iAx,jAt).  The  random  choice 
approximations  u  =  u(x,t;Ax,At)  satisfy  the  exact  equations  (1.1)  within  an 


error  of  the  form 


(7.7) 


£  E  (  i ,  j  )  =  //  4>  u  +  <|>  f  (u)dxdt 
(i, j)  X 


where 


(7.8) 


E(  i,  j )  =  /  $(x,(j  +  1 )  At)  {u  -  R(x/At,-u  ,u  )}dx 

.  n  we 

-Ax 

u  =  R(y  ...Ax/Atju  ,u  ) 
n  j+1  w  e 


It  has  been  shown  by  Liu  [26]  that  for  each  equidistributed  sequence  the  sum  of 
all  the  errors  E(i,j)  associated  with  mesh  diamonds  H(i,j)  vanishes  for  each 
given  test  function  <p  as  the  mesh  length  approaches  zero.  The  additional 
deviation  produced  by  the  above  relocation  of  certain  elements  yk  is  on  the 
order  Ax( |c^ |  +  !E2!)  an^  satisfies 

Ax( | |  +  | e2 I )  <  AxL|e1e2l/eo 

It  follows  that  the  sum  of  all  such  terms  vanishes  as  Ax  approaches  zero 


provided  that  o(Ax)  satisfies 


(7.9) 


lim  Ax/o(Ax)  =  0 
Ax+0 


Here  we  have  appealed  to  the  fact  that  the  total  amount  of  wave  interaction  as 
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measured  by  D  in  a  given  random  choice  approximation  is  bounded  uniformly  with 
respect  to  the  mesh  length  (12,  13] ,  which  in  this  context  implies 

l  l£1£2I  *  const.  , 

where  the  summation  is  taken  over  all  diamonds  associated  with  a  relocated 
element  y^. 

In  treating  systems  of  two  equations  with  eigenvalues  of  the  form 
X2  =  -X  ,  we  shall,  in  the  remaining  sections,  refer  to  the  random  choice 
generator  which  employs  an  equidistributed  sequence  satisfying  (7.4)  and  subject 
to  the  above  alteration  based  on  (7.6)  as  the  modified  random  choice  generator. 


8.  HYBRIDIZED  SCHEMES. 

In  this  section  we  shall  describe  a  class  of  hybridized  methods  based  on 
the  tracking  of  waves  whose  magnitude  lies  between  specified  thresholds 
depending  on  the  mesh  length.  Let  us  begin  by  considering  two  generating 
functions  based  on  the  same  grid  having  a  centered  diamond- shaped  stencil, 

(8.1)  u  =  <{>(u  ,u  ,u  ) 

n  w  s  e 

(8.2)  u  =  r(u  ,u  ) 

n  we 

where  <(>  corresponds  to  a  (possibly  perturbed)  scheme  in  class  N  and  r  to 
the  random  choice  method;  the  dependence  of  r  on  the  equidistributed  sequence 
is  suppressed.  Suppose  one  is  given  a  mesh  function  v  which  satisfies  either 
(8.1)  or  (8.2)  at  each  mesh  diamond,  i.e.  the  value  at  the  north  vertex  of  each 
mesh  diamond  is  obtained  from  the  three  lower  vertices  by  either  (8.1)  or  (8.2) 
according  to  some  presently  unspecified  rule  of  selection.  Let  C  and  RC 
denote  the  sets  of  diamonds  at  which  (8.1)  and  (8.2)  are  employed  and  consider 
the  associated  pattern  of  wave  magnitudes  in  v  as  introduced  in  Section  3. 

For  a  given  mesh  diamond  SI  let  SI  and  SI  denote  respectively  the  diamonds 
whose  ne  side  coincides  with  the  ws  side  of  Q  and  whose  wn  side 
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coincides  with  the  se  side  of  ft;  thus  an  incoming  wave  of  ft  crossing  its 
ws  side  serves  as  an  outgoing  wave  of  ft  while  an  incoming  wave  crossing 
its  se  side  serves  as  an  outgoing  wave  of  ft+.  Suppose  one  is  given  constants 
m  >  k  >  1  and  a  positive  function  a  =  o(Ax)  which  vanishes  in  the  limit  as 
Ax  approaches  zero.  In  terms  of  these  parameters  we  define  a  major  j-wave 
in  v  to  be  an  order  tuple 


/  1  2 

'VV 


whose  elements  e .  denote  the  magnitude  of  an  incoming  j-wave  of  some  diamond 
in  RC,  say  ft^  and  satisfy  the  following  conditions. 

(8.3)  If  el  is  the  outgoing  wave  of  a  diamond  in  C  then  |e.|  >  mo.  If 

1  3 

ej  is  the  outgoing  wave  of  a  diamond  in  RC  then  |e^|  >  ko  and  ft1 

contains  no  incoming  j-wave  6  such  that 

sgn  6_.  =  sgn  and  1 6_,  |  >  ko  . 
k  k+1 

(8.4)  If  p  >  1  then  £_.  and  e  serves  as  the  incoming  and  outgoing 

j -waves  of  ft^  and  satisfy 

k  k+1  -  k . 

sgn  e.  =  sgn  e,  and  |e.|  >  o  . 

3  j  3 

We  shall  say  that  a  mesh  diamond  ft  contains  an  incoming  major  j-wave  if  any  of 
the  following  five  conditions  hold. 


1) 

'V 

> 

m  a 

and 

ft' 

e  d 

2)  lYjl 

>  mo 

and 

ft+  E  D 

+ 

3) 

'V 

> 

k  a 

and 

ft' 

€  RC 

4)  l-r.l 

>  ko 

and 

ft  E  RC 

5) 

Either 

or 

Y3 

lies 

on  a  major 

j-wave 

We  note  that  if  any  of  the  first  four  conditions  hold  then  it  is  possible  for 
ft  to  serve  as  the  initial  diamond  of  a  major  j-wave. 

In  this  paper  we  shall  establish  the  stability  and  convergence  of  certain 
hybridized  schemes  whose  switching  operators  are  based  on  the  tracking  of  major 
waves.  The  starting  values  on  the  first  two  time-levels  t  =  0  and  t  =  At 
for  such  schemes  can  be  obtained  by  either  setting 
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u(jAx,0)  =  Ug{jAx);  u(kAx,At)  =  u^(kAx) 
for  appropriate  indices  j  and  k  or  by  using  standard  procedures  to  produce 
more  accurate  values  at  the  second  level  t  =  At  from  the  initial  data 
Uq(x).  Given  a  selection  of  parameters  m  >  k  >  1  and  o  =  a (Ax)  >  0  together 
with  starting  values  on  the  first  two  levels,  the  difference  approximation  is 
marched  forward  by  using  (8.1)  if  fl  contains  an  incoming  major  j-wave  for  some 
index  j  =  1,2,... ,n  while  (8.2)  is  used  otherwise.  For  a  class  of  systems  of 
two  equations  and  initial  data  with  small  total  variation  we  shall  establish  in 
Section  11  convergence  of  certain  hybridized  schemes  of  the  above  type  which 
employ  perturbations  of  class  N  generators  for  the  function  $  in  (8.1). 

It  would  also  be  interesting  to  study  hybridized  methods  based  on  a  single- 
level  threshold  where  the  random  choice  generator  is  used  to  compute  the  value 
at  the  north  vertex  of  a  mesh  diamond  if  there  exists  an  incoming  wave  <5_.  or 
Y ^  such  that  |6^|  >  o  or  Iy^I  >  0,  while  a  conservative  generator  <t>  is 
used  otherwise.  Such  schemes  correspond  to  more  standard  generating  functions 
of  the  form 

e<t>  +  (1  -  0)r 

for  an  appropriate  switching  function  9.  The  reason  for  introducing  a  multi¬ 
level  threshold  in  this  paper  is  to  ensure  the  convergence  of  the  corresponding 
difference  approximations  to  an  exact  solution.  Each  difference  approximation 
produced  by  hybridizations  of  the  above  type  satisfies  the  exact  equations  (1.1) 
within  an  error  which  consists  of  three  parts:  the  first  is  associated  with 
diamonds  in  C  and  vanishes  by  the  usual  arguments  employed  for  standard 
schemes  [22] ,  the  second  is  associated  with  the  equidistributed  sequence  on 
diamonds  in  RC  and  vanishes  by  the  argument  of  Liu  [26] ,  and  the  third 
involves  the  total  magnitude  of  waves  crossing  the  boundary  between  C  and 
RC.  In  Section  11,  we  shall  show  that  the  third  contribution  vanishes  as  Ax 
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approaches  zero.  To  do  this,  it  seems  necessary  to  rule  out  rapid  switching 


between  the  two  generators  (8.1)  and  (8.2)  during  the  propagation  of  an 
individual  wave.  In  hybridizations  of  the  above  type,  the  random  choice 
generator  remains  locked  on  a  major  wave  until  a  measurable  amount  of 
interaction  on  the  order  of  a(Ax)  has  reduced  its  magnitude  below  a.  The 
expectation  that  such  methods  will  produce  only  a  finite  amount  of  "wave 
interaction"  similar  to  that  of  the  exact  solution  operator  leads  one  to  suspect 
that  the  boundary  contribution  will  vanish  in  the  limit.  The  estimates  which 
make  these  remarks  precise  are  presented  in  Section  11.  In  this  connection  we 
also  recall  from  Section  7  that  we  must  require  o(Az)  to  vanish  slower  than 
Ax  in  order  that  our  modification  of  the  equidistributed  sequence  on  selected 
random  choice  diamonds  induces  a  deviation  from  the  standard  random  choice  error 
which  vanishes  in  the  limit.  This  restriction  also  prevents  a  premature  switch 
from  (8.1)  to  (8.2)  in  the  computation  of  a  focusing  compression  wave:  if  o 
has  the  same  order  as  Ax  then  the  random  choice  generator  is  engage  at  a  time 
on  the  order  of  At  prior  to  the  time  of  focus. 

Finally,  we  remark  that  it  would  be  interesting  to  study  the  corresponding 
hybridized  methods  which  are  based  on  incoming  major  j-waves  which  are  of  shock 
type.  For  technical  reasons  we  have  enlarged  the  class  of  incoming  major  j-wave 
to  include  both  shocks  and  rarefaction  waves.  We  conjecture  however  that  as  the 
mesh  is  refined  the  random  choice  method  is  in  fact  only  engaged  for  a 
substantial  length  of  time  on  the  major  shock  waves. 

9.  HYBRIDIZED  FUNCTIONALS. 

In  this  section  we  shall  construct  potential  functinals  which  will  he  used 
to  estimate  the  total  norm  for  the  hybridized  schemes  of  Section  7.  To  this  end 
we  first  eliminate  from  Pj(X)  the  terms  associated  with  pairs  of  j-rarefaction 
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waves,  since  they  are  not  compatible  with  the  random  choice  potential  which 
involves  only  pairs  of  approaching  waves.  We  redefine  the  weights  $ Ja.u^)  of 


(5.3)  by  setting 


(9.1) 


>.(a,u  )  =  b.(u  )  or  <j>.(a,u  )  =  b.(u  )  -  1/2 
3  a  3  a  3  a  3a 


if  a  is  a  j-shock  wave  or  a  j-rarefaction  wave  respectively,  crossing  the 
ws  side  of  a  mesh  diamond  with  base  state  and  by  setting 


(9.2) 


!>.(a,u  )  =  a.(u  )  or  4>.(<*»u  )  =  a.(u  )  -  1/2 
j  o  3  a  3  a  3  a 


if  a  is  a  j-shock  wave  or  a  j-rarefaction  wave  respectively,  crossing  the 

se  side  of  a  mesh  diamond  with  base  state  u  .  We  then  put 

a 


(9.3) 


P  .00  =  l  !a0|  +  l  <f>  (a,u  )a  , 
A  J  a 

j 


where  Aj  =  A^(X)  denotes  the  set  of  all  pairs  of  approaching  j-waves  in  X. 

As  before  the  functions  a^  and  b^  are  chosen  to  satisfy  (5.4)  and  we  define 


P(X)  =  l  P.(X)  . 
j=1  3 


The  behavior  of  such  functionals  on  local  wave  interactions  can  be 
efficiently  described  in  terms  of  the  connected  polygonal  arcs  which  consist  of 
line  segments  joining  adjacent  mesh  points.  Following  the  standard  convention 
[12],  we  shall  refer  to  such  an  arc  as  an  I-curve  if  the  x-component  varies 
monotonically  and  write  j2  >  J1  if  J 2  lies  toward  larger  time.  Two  I-curves 
J2  >  J.|  are  called  consecutive  if  they  coincide  except  along  the  boundary  of  a 
mesh  diamond  0  :  consists  of  the  wn  and  ne  sides  of  SI  while 

J1  -  J2  consists  of  the  ws  and  se  sides  of  S).  Here  it  is  natural  to 
associate  with  an  I-curve  J  of  a  given  mesh  function  u,  the  global 


configuration 


X  =  X ( J)  =  {(6k,Yk,uk)} 
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which  employs  the  values  of  a  set  of  alternating  vertices  pk  of  J  together 

K  k 

with  the  corresponding  magnitudes  of  interpolated  waves  6  and  y  .  In  the 
case  of  consecutive  curves  the  alternating  vertices  of  J2 

satisfy 

q.  =  p.  if  j  <  k  -  1  or  if  j  >  k  +  1 
j  D 

qk  =  n 

where  n  denotes  the  north  vertex  of  the  separating  diamond  ft.  Figure  1  of 
Section  4  can  be  used  to  illustrate  two  consecutive  I-curves  >  J1  as 
follows:  J1  passes  through  the  symbols  afc_v  uk'  sk'  uk+1'  sk+1  whH-e  J2 
passes  through  the  symbols  s}£_1 ,  u^,  w^ ,  uk+1,  sk+1.  Here  J.,  and  J2  are 
separated  by  the  diamond  ft  whose  vertices  are  labelled  by  the  symbols  u^,  sk , 
uk+1'  wk*  The  alternating  vertices  of  J1  are  labelled  by  the  symbols  sk 
while  the  alternating  vertices  of  J2  contain  the  sequence  sk^2'  sk-1'  wk' 

sk+1'  sk-2  * 

It  follows  easily  from  the  results  of  Section  7  that  the  new  functional 
P  satisfies 

(9.4)  P(J2)  “  <  "  od2(v)  +  c  osc(J  )C(V)  +  c|v|3 

if  TVJ ^  <<  1  and  if  and  J2  are  two  consecutive  I-curves  separated  by  a 

mesh  diamond  with  incoming  waves  V  =  ),  v,  =  (5.,y  ).  Here  we  put 

1  n  j  1  3 

2  «  2  n 
<3  (v)  =  l  (5  -  0  (u  )y  )  and  C(v)  =  £  C(6.,y.)  . 

j-1  3  3  8  3  j=1  33 

As  before  the  letter  c  denotes  a  positive  constant  depending  only  on  the 
equations,  the  scheme  and  the  state  u  in  R°  in  the  neighborhood  of  which  all 
analysis  takes  place.  The  proof  of  (9.4)  consists  of  the  observation  that  the 
term  1/2  introduced  in  (9.1)  and  (9.2)  accounts  for  the  previous  pairs  of 
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j -rarefaction  waves,  while  pairs  with  opposite  signs  introduce  only  a  quantity 
on  the  order  of  the  oscillation  times  the  cancellation. 

The  weights  employed  in  the  transverse  potential  do  not  require 
modification,  so  we  have 

T(J2)  -  T(J1)  <  -  ct{ n)  +  cd2(v)  +  c | v | 3  , 
for  the  perturbations  of  class  N  schemes  introduced  in  Section  6,  it  follows 
that,  for  approrpriate  constants,  the  functional  F  =  TV  +  cP  +  cT  satisfies 
F(J2)  -  F(J^ )  <  -  c{ T ( v)  +  d2(v)  +  C( v)}  +  c | v | 3 

if  TVJ1  «  1. 

Appropriate  functionals  for  the  hybridized  schemes  introduced  in  Section  8 
can  be  constructed  by  switching  the  weights  and  \Ji  of  T  and  P  on  and 

off  according  to  the  following  rule.  Associate  with  an  interpolated  wave  a 
the  unique  mesh  diamond  ft  for  which  a  is  outgoing.  If  ft  e  C,  let 
4> .  (a,u  )  keep  its  original  value  as  defined  by  (9.1)  and  (9.2),  otherwise  put 
<j>_.  =0.  If  a  and  0  cross  the  same  side  of  ft  and  ft  e  C  let  \|i(a,0) 
keep  its  original  value,  cf.  Sections  4  and  5,  otherwise  put  ij>  =  0.  Thus,  the 
weights  and  ip  vanish  if  and  only  if  the  associated  waves  are  outgoing  for 

a  random  choice  diamond.  Without  confusion,  we  shall  employ  the  notation  Pj 
and  T  for  the  hybridized  functionals  with  aforementioned  truncated  weights. 

The  stability  proof  for  the  hybridized  functionals  given  in  Section  10  is 
motivated  in  part  by  several  facts  concerning  the  relative  sizes  of  the  weights 
associated  with  the  physical  and  numerical  potentials  for  schemes  in  class  N. 
These  facts,  which  we  shall  present  below,  are  of  independent  interest  and 
hopefully  will  be  useful  in  future  investigations.  Let  us  now  restrict 


attention  to  systems  of  two  equations  with  eigenvalues  X^  and  X^  satisfying 
(1.2)  and  the  symmetry  condition  X^(u)  =  -X2<u).  The  well-known  system 


I 


0 


wt  +  p(v)x  =  °'  vt  "  wx  = 
arising  in  fluid  dynamics  and  elasticity  provides  an  example  of  symmetric, 

genuinely  nonlinear  eigenvalues  under  the  conditions  p'  <  0,  p“  >  0.  More, 

generally  we  note  that  systems  of  the  form  (1.4)  have  symmetric  eigenvalues 

X^  =  -X^ •  For  such  systems  it  is  natural  to  use  class  N  schemes  which  are 

symmetric  in  the  sense  that  they  employ  a  centered  stencil,  =  -0^ 

and  functions  w^  =  w^ta)  which  satisfy  w1  =  w2>  cf.  Section  3.  We  recall 

that  the  Lax-Friedrichs  scheme  is  symmetric  with  w^  =  w^  =  0.  A  brief 

calculation  show  that  for  such  symmetric  systems  the  weights  r  =  r12  and 

p  =  Pi j  associated  with  transverse  group  interactions  coincide  and  exceed  the 

weight  associated  with  approaching  transverse  wave,  i.e. 

r  =  p  >  s  . 

We  recall  that  if  w  =  0,  r  and  p  are  uniquely  determined  by  s  and  the 
local  base  state  according  to  formulas  (5.9).  For  such  symmetric  systems  we 
also  obtain  a  simplified  formula  discribing  the  effect  of  a  local  interaction 
onthe  transverse  potential  T: 

t(j2)  -  t(j1)  <  -  s|v1iiv2i(01  -  e2 )/( 1  +  e 1 ) ( 1  +  e2) 

+  (s  -  2p)ty^2  +  c |  v |  3  , 

where  <p.  =  6.  -  9.(u  )y  . .  The  formula  (9.5)  is  useful  in  connection  with  the 

j  3  3  s  j 

problem  of  forming  a  linear  combination  of  P  and  T  which  is  decreasing  on 
the  complement  of  the  second  order  weakly  interacting  states  W2»  To  this  end 
one  might  proceed  by  trying  to  dominate  the  indefinite  term  of  T, 

(9.6)  (s  -  2p)  \J»1  i{>2  , 

with  the  corresponding  negative  definite  terms  by  which  the  functional  P  is 
reduced,  i.e.  the  terms  under  the  summation  sign  in  the  right  hand  side  of  the 
following  inequality. 
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(9.7) 


P(J2)  -  PfJ^  <  ~  l  (Pj  +  t.)(1  -  PjMbj  -  V2)^ 

.  3 

+  c(osc  J  )C(v)  +  c|v| 

The  formulas  (9.5)  and  (9.7)  suggest  that  one  determine  the  restrictions  on  the 
weights  s,p  and  of  symmetric  schemes  which  imply  that  the  following 

quadratic  form  is  positive  definite: 

V  2 

(9.8)  l  (p.  +  T.)(1  -  p.)(b  -  1/2)4;.  +  (s  -  2p)4>-4'2  • 

3  3  3  3  3  1 

Here,  all  of  the  coefficients  are  evaluated  at  the  corresponding  local  base 
state  and  the  values  of  4» ^  are  regarded  as  arbitrary.  Now  for  a  fixed  choice 
of  s,  the  form  (9.8)  is  clearly  non-negative  for  appropriate  choices  of  b j . 
However,  a  brief  calculation  show  that  a  positive  definite  form  can  not  be 
achieved  with  a  weights  satisfying 

(9.9)  p/2  >  a1  -  1/2,  p/2  >  (b2  -  1/2)  . 

This  fact  is  of  particular  interest  in  connection  with  hybridized  schemes  since 
the  quantities  a^  -  1/2  and  b2  -  1/2  represent  the  maximal  weights 
associated  with  the  numerical  self-interactions  of  rarefaction  waves  as 
registered  by  we  have 

al  >  b1  and  b2  >  31 

since  (5.4)  holds  and  01  <  1  <  d^.  As  we  shall  see  in  Section  10  the  numerical 
group  interactions  characterized  by  coefficients  p  and  r  generally  lead  to 
favorable  contributions  in  bounding  the  hybridized  functionals  and  it  is  natural 
to  inquire  into  the  extent  to  which  they  might  compensate  for  the  less  favorable 
effects  of  switching  on  and  off  the  coefficients  associated  with  self¬ 
interactions. 
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Lastly,  we  remark  that  for  the  purpose  of  dominating  the  indefinite  term 
(9.6)  of  T,  one  might  appeal  to  the  favorable  term  with  leading  coefficient 
-s  on  the  right  hand  side  of  (9.5)  and  ask  what  restrictions  on  the  weights 
s,p  and  bj  of  symmetric  schemes  are  implied  by  the  condition  that  the 
following  form  be  positive  definite: 

l  (P.  +  T.)(1  -  p.)(b.  -  1/2)^2 
3D  3  3  3 

(9.io)  +  S(e2  -  e1  -  e) | vt 1 1 v2 1/( i  +  +  e2)  +  (s  -  2?)^^  . 

Here  the  quantity  e  >  0  is  introduced  simply  to  guarantee  that  a  residual  term 
of  order  t(v)  is  available  to  compensate  for  the  corresponding  growth  in  the 
total  variation  norm.  A  brief  calculation  shows  that  for  that  Lax-Friedrichs 
scheme,  and  hence  for  symmetric  schemes  with  w^  =  w2  sufficiently  small,  a 
positive  definite  form  (9.10)  can  be  achieved  with  a  choice  of  weights 
satisfying  (9.9).  Thus,  for  class  N  schemes  close  to  the  Lax-Friedrichs 
schemes  one  can  work  with  functionals  such  that  the  maximal  weights  associated 
with  numerical  self-interactions  of  shocks  are  less  than  half  the  weight 
associated  with  transverse  group  interactions. 

10.  STABILITY  OF  HYBRIDIZED  SCHEMES. 

In  this  section  we  shall  restrict  our  attention  to  systems  of  two  equations 
with  eigenvalues  satisfying  (1.2)  and  A^u)  =  -A2(u).  Our  results  extend  to 
the  more  general  class  with  eigenvalues  satisfying  (1.2)  and  A^  <  0  <  A2  but 
we  shall  treat  just  the  symmetric  case  A^  =  -A  for  concreteness.  We  shall 
establish  stability  in  the  total  variation  norm  for  hybridized  schemes  which  are 
based  on  the  tracking  of  major  waves  defined  by  parameters  k,m  and  o(Ax), 
cf.  Section  8,  and  which  employ  the  modified  random  choice  generator,  cf. 

Section  7  together  with  a  generator  of  the  form  $  +  q  where  $  corresponds  to 


r 


a  symmetric  scheme  in  class  N,  cf.  Section  9,  with  w^  and  Wj  sufficiently 
small  and  where  q  is  the  associated  perturbation  constructed  in  Section  6. 
Theorem  10.1.  If  a(Ax)  satisfies  (7.6)  then  for  appropriate  choices  of  m 
and  k,  m  >  k,  the  hybridized  schemes  above  produce  difference  approximations 
u(x,t)  =  u(x,t;Ax,At)  with  the  following  property.  If  the  initial  data  u0 
lies  in  a  small  neighborhood  of  a  state  u  e  Rn  and  if  TVUq  <<  1  then 

(10.2)  TVu(  ,t)  <  const. TVUg 

00 

(10.3)  /  |u(  ,t1 )  -  u(  ,t^ ) I dx  <  const. ( |t^  -  t  |  +  Ax)  , 

"  00 

2 

provided  that  t,t.j  and  tj  are  less  than  const. Ax/a  .  Furthermore,  the 
constants  depend  only  on  u,  the  equations  and  the  parameters  which  define  the 
scheme . 

Remarks .  The  approximate  L1  Lipschitz  continuity  (10.3)  for  the  marching  map 

is  an  immediate  consequence  of  (10.1);  the  difference  approximation 

u  =  u(x,t;Ax,At)  is  constant  on  each  rhombus  centered  at  a  mesh  point 

2  2  1/2 

(iAx,jAt)  with  sides  of  length  As  =  (Ax  +  At  )  oriented  by  the  normals 

a  and  0  and  assumes  therein  the  value  of  the  corresponding  grid  function  at 

(iAx,jAt).  One  may,  of  course,  derive  piecewise  constant  or  piecewise  smooth 

difference  approximations  from  a  given  grid  function  in  any  of  several  standard 

ways  and  still  maintain  estimates  of  the  form  (10.2)  and  (10.3).  We  shall 

comment  further  on  this  point  in  Section  11.  Lastly,  we  note  that  the 

particular  choice  o  =  Ax^  E  leads  to  uniform  estimate  (10.2)  and  (10.3)  over 

1-2e 

time  intervals  of  length  1/Ax  ;  the  ratio  of  meshlength  is  held  fixed  and 
satisfies  the  C-F-L  condition. 

Before  proving  Theorem  (10.1)  several  preliminary  remarks  are  in  order 
concerning  the  local  action  of  the  hybridized  functionals.  For  notational 
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convenience  we  shall  work  with  a  functional  of  the  form 

F(X)  =  cTVX  +  P  +  T 

where  c  is  an  appropriate  small  constant  and  P  and  T  denote  potentials 
with  weights  hybridized  according  to  the  rules  of  Section  9.  In  order  to 
distinguish  various  local  changes  in  F  according  to  the  structure  of  the 
incoming  waves  we  shall  employ  the  notation 


6.  =  6. 

3  3 

if 

n”  e  RC; 

S~  = 

3 

0 

if 

n  +e  c 

Y  .  =  Y  . 

3  3 

if 

e  RC; 

^3  = 

0 

if 

n+  e  c 

6.  =  <5. 

3  3 

-  6. 
3 

;  7.  =  Y. 

-73 

9 

if  6^  and  y.  are  incoming  j-waves  of  a  mesh  diamond  SI  whose  neighbors  with 
a  common  boundary  along  the  ws  and  se  sides  are  denote  by  SJ  and  ft+ 
respectively.  Given  two  consecutive  I-curves  separated  by  a  mesh 

diamond  SJ  with  incoming  waves  v  =  (6,y)  and  local  base  state  ug,  let 
P( v)  denote  the  incoming  potential  for  self-interactions  and  G(v)  the 
incoming  group  potential: 

2 

P(v)  =  l  $.(6.,u  )5.  +  ij>.(y.,u  )y . 

3  3  s  3  J  1  s  j 

G(v)  =  p(u  )|5  6  |  +  r(u  )|y  y  |  . 

S  1  Z  Si  z 

We  note  that  r  =  p  since  we  are  dealing  with  symmetric  schemes.  Now,  if 
TVJ1  <<  1,  we  have  the  following  estimates  with  an  appropriate  small  constant 
e : 

(10.4)  F( )  -  F(J2)  <  -(1  -  e)D(v)  -  eC ( v )  -  p(v)  -  G(v) 

(10.5)  -  FlJ2)  <  ~  eT<v>  “  Ed2 ( V)  +  P(V)  +  G(v)  +  0( | V | 3 ) 

if  S)  lies  respectively  in  RC  for  (10.2)  and  in  C  for  (10.3)  where 
D(v)  =  s  1 62  y^  |  +  J  {|6^y_,|  :  and  approach}  . 

The  estimates  (10.4)  and  (10.5)  can  be  summarized  by  introducing  a  quantity 
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I  (ft),  which  equals  the  negative  of  the  right  hand  side  of  (10.4)  if  !!  e  SC 
and  the  negative  of  the  sum  of  the  first  four  terms  on  the  right  hand  side  of 
(10.5)  if  ft  e  C,  and  obtaining 

F(J2)  -  F^)  <  I(ft)  +  const. o2(Ax)  ( v ( 

by  appealing  to  the  fact  that  the  waves  in  conservative  diamonds  have  an  order 
at  most  a.  Thus,  the  problem  of  proving  that  F  has  at  most  mild  growth  along 
orbits  requires  an  estimate  on  the  changes  in  F  due  to  the  switches  between 
generating  functions.  If  ft, ft  and  ft+  all  beong  toeither  C  or  if  ft  e  RC 
then  I  (ft)  <  0.  Hence  the  only  unfavorable  contribution  arises  if  ft  lies  in 
the  set  K  of  mesh  diamonds  in  C  with  at  least  neighbor  ft  ,ft+  in  RC.  If 
ft  e  K  then  F  is  augmented  by  the  effect  of  switching  on  the  coefficients 
associated  with  group  and  self-interactions.  Me  shall  show  that  for  each 
diamond  ft  in  K  one  can  associate  a  set  of  diamonds  at  lower  time  levels  at 
which  compensating  interactions  have  occurred.  Specifically,  we  shall  show  that 
if  TVMkX  is  sufficiently  small  for  k  <  p  then 
(10.7)  l  {i (ft)  :  ft  [0,kAtJ }  <  0  , 

if  k  <  p.  Combining  (10.5)  with  the  property  that 

F(MX )  <  F(X){1  +  const. 02(Ax)}  -  £  I(ft.) 

j  3 

where  {ft_.}  denotes  the  set  of  all  mesh  diamonds  separating  the  I-curves 
associated  with  X  and  MX,  we  obtain 

F(M^X)  <  2F(X)  if  p  <  const. /o2 (Ax) 

and  TVX  <<  1.  Furthermore,  in  the  course  of  the  proof  of  (10.7)  we  shall 
establish  a  slightly  stronger  estimate  which  facilitates  the  proof  that  the 
difference  approximations  converge  to  an  exact  solution.  The  stronger  estimate 
assumes  the  form 
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( 10.8) 


l  1(B)  <  -  £  I  I*(ft) 


where  0  <  e  <<  1,  both  summations  are  taken  over  the  mesh  diamonds  ft  which 
lie  in  a  given  strip  [0,pAt]  and  the  interaction  term  1*  is  defined  by 

(10.9)  I*(ft)  =  t(v)  +  d2(v) 

(10.10)  I*(ft)  =  D( V )  +  C(v)  +  G(v)  +  PQ(v)  +  H ( v ) 

if  ft  lies  respectively  in  C  for  (10.9)  and  RC  for  (10.10).  Here  PQ 
records  the  self-interactions  associated  with  those  incoming  waves 
v  =  (6/Y)  of  ft  which  have  maximal  coefficients  and  enter  from  a  diamond  in 


C,  i . e . 


V”>  ■  *  t2  ■ 

and  H(v)  records  the  effect  of  switching  off  coefficients  at  the  starting 
points  of  major  waves: 


2 

H ( v)  «  mo 

if  ft  contains  an  initial  segment  of  a  major  wave  while  H=0 ,  otherwise.  The 
factor  of  m  is  introduced  merely  as  a  reminder  of  the  definition  of  H.  We 
deduce  from  (10.8)  that 

F(M^X)  <  2F(X)  -  const.  £  I*(ft)  if  k  <  const. /o2(Ax) 
and  therefore  that  total  amount  of  interaction  as  measured  by  I*  in  a  strip  of 
the  form  [0,pAt]  satisfies 


(10.11)  £  {l*(ft)  :  ft  e  [0,pAt]}  <  const. TVuQ 

2 

if  p  <  const. /o  . 

Proof  of  Theorem.  Consider  a  fixed  strip  [0,pAt]  and  denote  by  K  the  set  of 
mesh  diamonds  ft  therein  such  that  either  ft  or  ft+  lies  in  RC.  We  shall 
indicate  how  to  associate  with  each  ft  in  K  a  collection  of  diamonds 
r  =  T^(ft),  j  =  1,2, ...,£,  and  compensating  quantities  q(I\,ft)  such  that 
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i  -  *  .  , 


It  then  follows  that 


l  usu  +  l  i (ft)  <  l  mn  +  l  l  q(f.,n)  <  o  . 

K  Kc  K  K  j  3 


Here  diamonds  of  the  form  T.(Q)  and  T  (ft')  may  coincide  for  distinct 

1  k 

ft  and  ft*.  The  estimate  (10.6)  simply  follows  from  the  fact  that  compensating 
quantities  q  can  be  located  without  appealing  to  a  small  fraction  of  I*(ft). 

Fix  SI  in  K.  In  the  case  study  below  we  shall  locate  diamonds  I\ 
having  interactions  which  compensate  for  the  contribution  to  I (ft)  from  waving 
leaving  T+.  The  set  of  all  diamonds  of  the  form  I\  or  provides  the 

desired  collection  {I\},  Thus,  let  us  suppose  that  ft  lies  in  RC  and  let 
<S  and  S  be  the  waves  in  the  Riemann  problem  with  data  (u  ,u  )  where  w 

1  4b  W  6 

and  e  denote  the  values  of  the  difference  approximation  at  the  west  and  east 
vertices  of  ft  .  We  shall  first  treat  the  case  of  two  shock  waves.  The 
remaining  cases  are  somewhat  easier  and  will  be  discussed  below. 

Case  1 .  <S  ^  <  0,  62  >  0 

Subcase  1 .  Suppose  both  6^  and  ^  cross  the  ne  side  of  ft  .  If  ft 
contains  an  incoming  major  j-rarefaction  wave  then 

C(ft")  +  D(ft“)  >  CO 

since  there  exists  no  outgoing  rarefaction  wave  for  ft  .  As  before  we  shall 
denote  by  c  any  of  various  constants  depending  only  on  the  scheme,  the 
equations  and  the  neighbor  in  which  all  analysis  takes  place.  On  the  other  hand 


I( ft)  <  ck^o^  . 

Hence,  in  this  situation  we  associate  with  ft  just  a  single  diamond  r  =  ft 
and  put 

q(r  ,ft)  =  OC(ft")  +  CTD(ft~)  . 

Now,  if  ft  contains  an  incoming  major  j-shock  then 

l  c(ftA)  +  D(fti)  >  ca 

where  {ft^}  denotes  the  set  of  diamonds  which  it  intersects.  Here,  we  put 
ri(ft)  =  ft^  and 

q(r±,ft)  «  0C(fii)  +  oDjft^  . 

Subcase  2.  Suppose  5^  crosses  the  wn  side  of  ft  while  62  crosses  the 

ne  side  of  ft  .  If  ft  contains  a  major  2-wave  the  analysis  of  subcase  1  is 

applicable.  Let  us  therefore  consider  the  situation  where  ft  contains  a  major 

1-wave.  Here  the  only  non-negative  contribution  to  I(ft)  associated  with  waves 
-  2 

leaving  ft  is  given  by  b2^2  an<3  aPProPriate  compensation  can  be  obtained  as 
follows.  Let  6  be  an  outgoing  j-shock  of  an  arbitrary  diamond  in  RC.  We 
introduce  a  backward  shock  tree  T(6)  through  6  by  starting  with  6  and 
repeating  the  following  process:  given  an  outgoing  j-shock  of  a  diamond  in 
RC  associate  with  it  the  corresponding  set  of  incoming  j-shocks,  if  any.  Such 
a  tree  T(6)  is  contained  in  a  union  of  RC  diamonds,  say  ft^,  each  having 
one  outgoing  j-shock  and  at  most  two  incoming  j-shocks.  With  reference  to  the 
particular  shock  $2  we  shall  refer  to  a  corresponding-  diamond  ft^  as  a 
terminal  diamond  for  T(62>  if  any  of  the  following  conditions  hold: 

(10.13)  ft^  lies  in  C  and  contains  an  outgoing  2-shock, 
say  c»2 ,  crossing  its  wn  side 

(10.14)  ftj^  lies  in  C  and  contains  an  outgoing  2-shock, 
say  S2  crossinq  its  ne  side 

(10.15)  ft,  contains  no  incoming  2-shock. 

k 
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We  shall  refer  to  such  associated  waves  a and  as  terminal  shocks  for 

T(62).  We  can  estimate  <S2  in  terms  of  terminal  shocks  of  T($2)  as  follows: 

(10.16)  <  I  a2  +  l  e2  +  2  ^  W  +  C  °SC  T(62)  E  d(V  * 

k  k 

Here  the  first  two  summations  are  taken  over  all  terminal  shocks  of  T(62>  as 

indicated  in(10.13)  and  (10.14);  S^ft^)  records  the  shock-shock  interaction  of 

the  second  field  in  ft,  ,  i.e. 

k 

S2(ftk)  *  | Eh | 

if  ft^  contains  incoming  2-shocks  £  and  n  while  S2  =  0  otherwise;  and 

osc  T(S  )  denotes  the  largest  wave  magnitude  in  the  union  of  ft  •  Thus,  in 

the  case  of  difference  approximations  with  small  oscillation,  the  sum  of  squares 

of  terminal  shocks  plus  the  total  shock-shock  interaction  on  T(6  )  bounds  the 
2 

"output"  62  modulo  a  small  fraction  of  the  total  random  choice  interaction 
occurring  on  T(<$2).  The  estimate  (10.16)  can  be  established  by  proving  by 
induction  that 

S2  ‘  l  «2  *  I  *  2  l  W 

+  4c  osc  T($  )  T  D(fi  )  +  c  y  D  (ft,  ) 

2  L  k  u  k 

where  the  constant  c  is  chosen  so  large  that  the  incoming  and  outgoing  j-waves 
of  a  general  random  choice  interaction  satisfy 

'Gj  ~  (Sj  +  <  cD( 5  ,y)  , 

cf.  Section  7.  At  this  point  we  remark  that  we  can  choose  the  weights  b2 
(and  a.j)  to  be  constant  since  we  only  need  (5.4)  to  hold.  If  b2  is 
constant,  we  obtain  an  estimate  of  the  form 

V2  ‘  £  b2a2  *  ‘  b2C2  *  =b2°SC  T<V  J  D<V 
*  2b2  J  W  ' 


(10. 17) 


for  which  there  exist  obvious  compensating  quantities  q(r,,f2)  for  the  first 

three  term  on  the  right  hand  side.  It  remains  only  to  obtain  compensation  for 

the  last  term  of  the  right  hand  side  of  (10.17).  Now,  since  bj  >  1/2,  the 

shock-shock  interaction  (with  coefficient  equal  to  one)  in  D((2,  )  does  not 

k 

cover  this  last  term.  However,  we  can  obtain  adequate  compensation  for  it  by 

appealing  to  diamonds  of  the  form  (2^.  Specifically,  in  subcase  2  we  shall 

associate  with  a  collection  of  diamonds  T .  such  that 

1 

ur.cua  u  n"  . 

j  k  k 

We  partition  the  analysis  of  subcase  2  as  follows. 

Subcase  2A.  Suppose  satisfies  (10.13).  Let  0^,0^  c*eno'te  the  outgoing 

waves  of  (2^  crossing  the  ws  side  of  (2^  and  ,a ^  the  outgoing  waves  of 
(2*  crossing  the  se  side  of  ft  .  Now  if  (2  contains  a  major  2-wave  then  we 

K  K  K 

may  proceed  exactly  as  in  subcase  1.  Let  us  assume  therefore  that  ft^  contains 

a  major  1-wave.  If  lies  on  a  major  wave  then  lo^l  >  ma  and  we  can  appeal 

2 

to  (a  small  fraction  of)  the  quantity  -a^a^  which  appears  in  Ifft^)?  we  have 

2  2  2 
~ai 01 1  <  -a,m  0  ' 

with  m  large.  Hence,  the  only  difficulty  arises  in  the  situation  where  we 
know  only  that  lies  on  a  major  1 -wave  and  this  can  be  handled  as  follows. 

Since  01  and  0^  cross  the  same  line  segment,  the  definition  of  the  modified 
random  choice  generator,  cf.  Section  7,  implies  that  either 

I  I  <  eo2  or  |  61  I  +  1 62 1  >  La  . 

2 

Suppose  the  former  inequality  holds.  We  write  the  term  a^ot^  appearing  in 

-I  (ft,  )  in  the  form 
k 

MO-10)  a2a2  =  b2a2  +  (a2  ‘  b2)0t2  * 

We  recall  that  aj  depends  on  the  value  of  the  difference  approximation  at  the 
south  vertex  of  ft^  and  satisfies  a^  >  b_>.  The  first  term  on  the  right  hand 
side  of  (10. 10)  appears  in  (10.17)  and  we  have  a  residual  of  the  form 


(«2  “  *52^a2  w*1*c*1  we  take  advantage  of  as  follows.  We  note  that  if 


(10. 19) 


(*2  “  b2,a2  >  2b2S2(V 


we  have  located  adequate  compensation  for  the  shock-shock  interaction  in 


Now  (10.19)  holds  if 


IS2I  <  C.2  -  b2)la2l/21>2 


since  =  ^2°2^*  Equivalently »  (10.19)  holds  if 


(10.20) 


CO  <  (a2  -  b2)|a2l/2b2 


since  (3^  lies  on  a  major  wave  and  <  eo  .  Therefore  no  further 

analysis  is  required  unless  we  are  faced  with  the  opposite  inequality  from 


(10.20),  i.e.  unless 


(10.21) 


l<*2 1  <  2b2eg/(a2  -  b2)  . 


However  if  (10.21)  holds  we  can  locate  adequate  compensation  q  in  the  diamond 
as  follows.  Inequality  (10.21)  implies  that 

VV  ‘  l2V/(*2  -  V1o62  • 

In  a  lemma  below  we  shall  show  that  |  0 ^  62  I  *  cl(fl  )  where 

i(n~)  =  g(«~)  +  D(n~)  +  H(n“)  . 

We  therefore  obtain  a  1 02  |  <  cKfi^)  and 


(10.22) 


2b2s2<V  <  Eci(n  )  . 


We  conclude  from  (10.22)  that  if  |  (3  8  |  <  E0  there  exists  an  appropriate 
compensating  quantity  q  for  the  shock-shock  interaction  in  since  we  have 


an  estimate  of  the  form 


(10.23) 


b2°2  +  2b2S2(V  <  EC?(V  +  a2a2 


with  small  e.  We  note  that  (10.13)  can  be  strengthened  by  the  inclusion  on  the 

2 

left  hand  side  of  a  term  of  the  form  This  can  easily  be  seen  by  writing 

(10.18)  in  the  form 

a2a2  “  b2a2  +  (a2  ~  b2,a2'/2  +  (a2  ~  b2)a2//2 
and  proceeding  as  above.  This  accounts  for  the  presence  of  the  term  Pn(v)  in 
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(10.10).  We  complete  the  discussion  of  subcase  2A  by  considering  the 
alternative  situation  where  18^  +  |  82  I  >  La.  Here,  if  |82I  >  ma  we  can 
appeal  to  the  cancellation  process  exactly  as  in  subcase  1.  On  the  other  hand 
if  | 82 I  <  mo  we  obtain 

1(0  >  c (L  -  m)o2 
k 

since 

1 61  e2  I  <  cl(£2~)  and  18^  >  (L  -  m)o 

Thus  if  L  is  sufficiently  large  it  follows  that  I(£2)  is  bounded  by  a  small 
fraction  of  I(£2^)  and  the  desired  compensation  is  realized. 

Next  we  shall  establish  the  lemma  required  in  subcase  2A. 

Lemma.  If  ei  and  denote  a  pair  of  outgoing  waves  crossing  either  the 

wn  or  ne  side  of  a  diamond  (2  in  RC  then 

f e1 e2 f  <  cl(£2)  5  c{G(£2)  +  D(£2)  + 

Proof :  Let  6.|,62  d®note  the  incoming  waves  crossing  the  ws  side  of  £2  and 
Y-jfYj  tfie  incoming  waves  crossing  the  se  side  of  £2.  We  have 

*Ej  ”  +  Yj) I  <  cD(Q)  . 

If  £2+  and  £2  both  lie  in  RC  the  structure  of  the  modified  random  choice 

generator  implies  that  51  =  =  0  and  hence 

*  E 1  e2  I  <  C^2Y1  *  =  * 

If  £2  e  C  and  £2+  e  RC  then  Y2  =  0  and  we  obtain 

|e1e2l  <  1 6<1 6  |  +  cD  (£2)  <  cG(£2)  +  cD(£2)  . 

Similarly  if  £2+  e  c  and  £2  e  RC  then 

|  e2 1  <  cG(£2)  +  cD(£2)  . 

Finally,  if  £2  and  £2  both  lie  in  C  we  obtain  the  simple  estimate 

2  2 

I e ^ e2 1  <  cm  a  <  cH(£2)  . 

The  proof  of  the  lemma  is  complete. 


Subcase  2B.  Suppose  that  ft^  satisfied  (10.14),  as  in  subcase  2A  the  situation 

if  simple  unless  we  know  only  that  ft^  contains  a  major  1-wave.  Supposing  this 

to  be  the  case,  we  proceed  as  follows.  If  8,  lies  on  a  major  1-wave  then 

2 

I8.J  >  mo  and  we  can  appeal  to  (a  small  fraction  of)  the  quantity  -h,^ 
which  appears  in  I(ft^)  and  obtain 

2  2  2 
-b^  <  -b^na 

for  appropriate  compensation.  Therefore,  let  us  assume,  on  the  other  hand,  that 

lies  on  a  major  1-wave  and  let  us  estimate  S2  ^  *  Now  S2^k^  vanishes 

if  a2  =  0.  If  a2  *  0  then  either 

2 

1  °ti  ot2  I  <  eo  or  lai  I  +  la2  I  >  La  * 

If  the  former  inequality  holds  then 

s2(V  =  I ct2 B2 1  <  eo|82l  <  elct^l  <  eD(nk)  • 
since  |a  |  >  a.  If  the  latter  inequality  holds  then  we  remark  that  either 
| a2 |  >  ma  in  which  case  we  can  appeal  to  the  cancellation  process  exactly  as 
in  subcase  1  or  | a2 |  <  mo  in  which  case 

Dink)  >  *e2ai*  ^  *02*<L  ~  m>°  >  (L  "  m>  f e2 a2  * /”* 

since  | ot ^  |  >  (L  -  m)o.  Hence  if  L  is  sufficiently  large 

“sVV  ‘ ;  D'V 

and  we  conclude  that  a  small  fraction  of  the  random  choice  interaction  in  ftfc 
provides  adequate  compensation  for  S^ft^)  (if  m  is  large). 

Subcase  2c.  Suppose  ft,  is  not  a  terminal  diamond.  Then  ft*  and  ft,  both 

—  1C  K  Jc 

lie  in  RC  and  we  proceed  as  follows.  If  8^  *  0  then  <*2  =  0  and  no  further 
analysis  is  necessary.  If  81  =  0  then  we  may  assume  that  Ic^l  >  o.  Now, 
if  a2  <  0  then  (10.24)  holds.  Under  the  first  inequality  of  (10.24)  we  obtain 
an  adequate  estimate  of  the  form 

JWV  ■  2b2la282'  ‘  eoD«V  • 


Suppose  the  second  inequality  of  (10.24)  holds.  If  | |  >  mo  then  no  further 
analysis  is  required  as  our  previous  cases  apply.  If  la^l  <  mo  then 
I a^|  >  (L  -  m)o  and  we  obtain 

D(&k)  >  slB^I  >  s(Ii  -  m)a|B2l  >  s(L  -  m)|a232l/2 
=  s(L  -  m){2b2S2(f2k)}/4b2  , 

since  we  may  assume  that  J ct2 1  <  2a,  otherwise  the  situation  if  trivial. 

Hence,  if  L  is  sufficiently  large,  a  small  fraction  of  0(5)^)  dominates  the 
shock-shock  interaction  of  the  second  field  in  This  completes  the  proof  of 

case  1 . 

We  note  that  the  presence  of  the  terms  (10.9)  and  (10.10)  in  I* (ft)  is 
simply  due  tothe  fact  that  we  were  able  to  obtain  appropriate  compensation  for 
the  positive  terms  in  I(fl)  by  appealing  to  only  a  fraction  (less  than  one)  of 
the  available  negative  terms. 

It  remains  to  discuss  the  cases  where  the  Riemann  problem  with  data 
(uw,ue)  produces  a  shock  and  a  rarefaction  wave  or  a  pair  of  rarefaction 
waves.  We  note  that  if  the  strip  under  consideration  involves  no  splitting  of 
rarefaction  waves  by  the  sample  point  then  the  analysis  of  case  1  applies 
virtually  without  modification  to  the  remaining  cases.  If  a  sample  point  splits 
a  rarefaction  wave,  only  a  small  change  in  the  above  analysis  is  required  since 
the  splitting  of  rarefaction  waves  is  not  by  itself  produce  any  wave 
amplification.  In  this  connection  we  note  that  it  is  primarily  to  handle  one 
technical  point  connected  with  wave  splitting  that  we  restrict  attention  to 


symmetric  schemes  with  w1  and  w2  small. 

We  begin  with  a  remark  needed  to  construct  the  analogous  rarefaction 


trees.  Let  and  denote  the  j-waves  crossing  the  wn  and  ne  side  of 


Case  C.  If  6^  >  0,  *  0  then  partitions 


s)  ‘  j,  V  “a 


can  be  formed  in  such  a  way  that  either 


(10.26) 


f  l*,k  -  <*,J  +  f  -  r.  .  I  <  cD(Si) ,  or 

k=1  3k  ]k  p+1  jk  3  »k-p 


no.27)  1  ^  ^,1  *  ",p-« 1  ‘  “«»>  • 


We  note  that  virtually  the  same  magnitude-preserving  partitions  are  employed  in 
[13]  for  the  purpose  of  constructing  characteristic  curves.  Using  such 
partitions,  an  analogous  backward  rarefaction  tree  through  a  given  outgoing 
rarefaction  wave  of  a  diamond  ft  in  RC  can  be  constructed  by  repeating  the 
following  process  which  associates  with  a  given  outgoing  j-rarefaction  wave 

tt  the  corresponding  incoming  wave  or  waves  according  to  the  following  rules 
j  * 


In 

Case 

A: 

"jk 

~  <5- 

3k 

In 

Case 

B: 

"jk 

«*  v  • 

jk 

In 

Case 

C: 

If 

(10.26)  holds  then 

V  ~  V  if  k  <  p  and  wjk  ~  Yj,k-P  if  k  >  p  +  1  • 

If  (10.27)  holds  then 

"jk  *  V  lf  "  ‘ p  - 1  *"a  ’jk "  Vk-p*.  if  * » » *  ’ 

’jk  *  1  VYj.>  • 

Thus,  only  in  Case  C,  (10.27)  does  not  associate  a  pair  of  waves  with  a  given 


outgoing  wave.  We  now  proceed  to  sketch  the  remaining  cases. 


crosses 


Subcase  1 .  Suppose  that  6^  crosses  the  vm  side  of  ft  while  <5^ 
the  ne  side  of  ft  .  This  subcase  can  be  handled  exactly  as  in  Case  1. 

Subcase  2 .  Suppose  that  6^  is  split  by  the  sample  point  into  a  1-wave  6^ 
crossing  the  wn  side  of  ft  and  a  1-wave  6”  crossing  the  ne  side  of  ft 

along  with  <5^ :  6,|=6'+6y.  Here  we  must  locate  compensation  for  terms  of 

the  form 

v;2  ♦  »s;s2 

which  enter  X(ft).  Now  if  ft  contains  either  an  incoming  major  2-wave  or  an 

incoming  major  1-shock  the  analysis  of  Case  1  suffices.  Let  us  therefore  assume 

that  ft  contains  an  incoming  major  1-rarefaction  wave.  Let  denote  the 

incoming  waves  of  ft  crossing  its  ws  side  and  a^cx^  *'*le  incomin9  waves  of 

ft  crossing  its  se  side.  Let  ft^  and  ft2  denote  the  diamonds  which 

intersect  the  ws  and  se  sides  of  ft  respectively.  For  concreteness  we 

treat  the  situation  where  ft^  and  ft^  both  lie  in  RC;  this  case  provides  the 

primary  example  of  the  switching-on  of  weights  due  to  the  splitting  of  a 

rarefaction  wave.  Here  we  proceed  by  tracing  back  the  major  2-rarefaction 

wave.  Here  we  proceed  by  tracing  back  the  major  2-rarefaction  wave  through 

to  the  diamond  T  which  contains  its  initial  segment  e.  Now,  if  e  enters 

T  from  a  diamond  in  C  then  e  >  mo  and  we  can  appeal  to  a  quantity  on  the 

2 

order  of  e  appearing  in  I(D.  On  the  other  hand,  if  e  enters  T  from  a 

diamond  in  RC  then  e  >  ka  and  we  shall  show  that  there  again  exists 

2 

compensating  interactions  on  the  order  of  e  in  certain  diamonds  near  F  to 
be  described  below. 

Before  proceeding  with  the  latter  case  we  shall  remark  on  the  situation 
where  more  than  one  major  rarefaction  wave  is  traced  back  to  the  same  diamond 
T.  Now,  if  {e^}  denotes  the  collection  of  all  waves  of  the  above  type  which 


can  be  traced  back  to  the  same  e,  then,  modulo  a  quantity  equal  to  a  small 


fraction  of  the  random  choice  interaction  in  the  associated  diamonds,  we  have 

e2  >  m  J  e2  if  e  >  mo  and  e2  >  k  Y  e2  if  e  >  ka  , 

L  j  L  j 

since  in  both  situations  <  a.  Indeed,  without  loss  of  generality  we  have 

E  “  l  ej 

modulo  interactions,  and  therefore  in  the  case,  say,  e  >  mo, 

e2  -  I  e2  +  e  J  e ^ d  -  e^/e)  >  l  e2  +  (1  -  1/m)e2  . 

Thus,  the  splitting  of  rarefaction  waves  produces  a  situation  where  the  sum  of 
the  squares  of  the  terminal  waves  associated  with  the  switching-on  of 

weights  is  a  small  fraction,  here  1/m  or  1/k  with  m  and  k  large,  of  the 
square  of  the  initial  wave  e,  modulo  a  quantity  on  the  order  of  the  same  small 
fraction  of  the  total  random  choice  interaction  as  measured  by  D  in  the  union 
of  diamonds  through  which  the  traced  waves  pass. 

Our  analysis  of  the  case  where  e  enters  T  from  a  diamond  in  C  is 
subdivided  as  follows. 

Subcase  2A .  Suppose  T  and  r+  both  lie  in  RC  let  (6^, ^5  and  (Y^Y^) 
denote  the  incoming  waves  of  T  crossing  its  ws  and  se  sides.  We  have 


(10.28) 


e  -  61  +  yi  +  0D(D  . 


Now  if  6^  »  Y1  =  0  then  |e|  <  cD(T)  and  we  easily  obtain  appropriate 

2 

compensation  of  the  form  e  by  appealing  to  a  small  fraction  of  D(F).  Indeed 


we  have 


e  <  c(osc  u)D(n  , 


since  D(T)  is  itself  on  the  order  of  the  oscillation  of  the  difference 


approximation  u.  Let  us  therefore  suppose  that  5^  and  y1  do  not  both 
vanish.  For  concreteness  assume  6^  *  0.  It  then  follows  from  the  definition 
of  the  modified  random  choice  generator  that  y2  “0*  If  Y1  =  0  then 


■  6^  contradicting  the  fact  that  e  is  the  initial  segment  of  a  major 
rarefaction  wave.  On  the  other  hand  if  y1  *  0  then  a  simple  calculation  shows 
that 

(10.29)  16^1  +  |Yll  <  c|62I  , 

since  both  61  and  y  are  non-zero  rarefaction  waves.  In  this  connection  we 

note  if  *  0  then  at  least  one  of  the  waves  61#Y.|  must  vanish  since  the 

associated  eigenvalue  increases  from  left  to  right  across  1-rarefaction 

waves.  It  follows  from  (10.28)  that 

D(D  +  l(r")  >  c|62y1I  +  c| 61 62 1  >  0(16^  +  IyJ)2 
and  we  conclude  that 

e2  <  c{d(D  +  l(r")}  , 

to  which  we  appeal  for  appropriate  compensation.  The  situation  where  y.j  *  0 
is  handled  in  a  similar  fashion. 

Subcase  2B.  Suppose  T  and  T+  both  lie  in  C.  If  (6^1  <  lY^f  we  appeal 
2 

to  -a^.,  appearing  in  I(T)  and  obtain 

e2  <  8Y2  +  2cD(T“) 
using  (10.28).  Hence  a  quantity  on  the  order  of 

-alY2  +  (osc  u)D(T) 

2 

bounds  e  .  If  16^  >  Iy^  we  proceed  as  follows.  Either  y2  or  &2  lies 
on  a  major  2-wave.  If  Y2  does  then 

(10.30)  t  <  2ka  +  cD<D 

2 

since  <  kc.  In  this  situation  we  appeal  to  the  quantity  “a2^2  aPPe*ring 

in  I ( r )  and  note  that  modulo  a  small  fraction  of  D(T)  we  have  an  estimate  of 


the  form 

2  2  2  2.2.  2 
*2^2  *  a2m  8  >  a2m  ^  E 

using  lYjl  >1*0.  Here  we  are  led  to  take  m  much  larger  than  k.  Next, 


suppose  that  62  lies  on  a  major  2-wave.  In  this  situation  we  simply  repeat 
the  argument  given  in  subcases  2A  and  2B  and  in  subcase  2C  below  with  T 
replacing  T  and  use  the  fact  that,  modulo  the  interaction  term  D(T)  we  have 
>  ka/2  since  16^1  <  Iy^. 

Subcase  2C.  Suppose  r  lies  in  C  while  T+  lies  in  RC.  If  |6.|  >  Iy.jI 
2 

we  appeal  to  appearing  in  I(T)  and  obtain  an  inequality  of  the  form 

b^2  >  b  e2/4  -  cD2(T)  . 

Next,  let  us  suppose  that  16^1  <  Iy^.  If  <$2  lies  on  a  major  wave  then  we 
2 

appeal  to  “b2^2  appearing  in  I(T)  as  follows.  We  have  (10.30)  since 

5^  <  ka  and  16^  <  Iy^.  Therefore 

o  2  2  2  2  2  2 

b262  >  b2m  0  >  b2m C /4k  ~  CD  (r)  * 

On  the  other  hand  if  y2  lies  on  a  major  2-wave  we  simply  repeat  the  arguments 
above  replacing  T  with  T  and  using  the  fact  that 

Y1  >  ko/2  -  cD(D 
This  completes  the  sketch  for  the  three  subcases. 

The  remaining  term  p6"62  can  be  h*ndled  as  follows.  From  the  definition 
of  the  modified  random  choice  generator  we  have  ■  82  =  0  since  ^  and 
lie  in  RC.  Thus,  a  term  of  the  form 

sla^l  +  ck6^'2 

is  available  to  compensate  for  p6^62.  Writin<?  “  TB1  with  0  <  t  <  1,  we 
require  an  inequality  of  the  form 

(10.31)  pl<S"«2f  =  PT  la2^i  I  <  s|0ia2l  +  ^  • 

But  since  t  S ^ /a2  I  *  Vk,  (10.31)  is  equivalent  to  the  condition 

2 

(10.32)  pT  <  s  +  ex 

Now  an  appropriate  restriction  on  the  C-F-L  number  guarantees  that  p  is  only 
slightly  larger  than  s  which  implies  that  (10.32)  is  valid.  We  lastly  remark 


that  the  case  where  6^  <  0  and  <5^  *  0  and  the  case  where  6^  >0  and 
62  >  0  are  treated  in  the  same  fashion  as  above. 

11 .  CONVERGENCE  OF  HYBRIDIZED  APPROXIMATIONS. 

In  this  section  we  shall  consider  the  class  of  systems  and  hybridized 
schemes  for  which  we  established  stability  in  Section  10  and  show,  under  the 
additional  hypothesis  that  the  error  function  i associated  with  the 
equidistributed  sequence  and  the  switching  function  o(Ax)  satisfy 

lim  i|»(n)n^^  =  0  and  lim  \|;(  1/Ax)/(  Axo)  ^  ^  -  0  , 
n+«  Ax+0 

that  a  subsequence  of  difference  approximations  converges  to  an  exact 
solution.  We  note  that  standard  arguments  using  the  stability  estimates  (10.2) 
and  (10.3)  yield  the  existence  of  a  subsequence  of  difference  approximations 
converging  pointwise  a.e.  to  a  function  u  which  is  a  function  of  bounded 
variation  in  the  sense  of  Tonelli-Cesari  [10,  33];  indeed,  the  subsequence 
converges  to  u  in  L1  of  the  space  variable  for  each  fixed  time  t.  We  also 
note  in  passing  that  our  stability  analysis  shows  that  the  total  variation  of 
the  difference  approximations  along  (space-like)  lines  with  speed  of  propagation 
greater  than  or  equal  to  Ax/At  in  absolute  value  is  uniformly  bounded  and 
consequently  that  there  exists  approximate  L1  Lipschitz  continuity  in 
directions  normal  to  such  lines.  It  remains  to  prove  that  u  satisfies  (1.1) 
in  the  sense  of  distributions,  equivalently  that  the  contour  inteqral 

/  v(u)ds;  v(u)  =  v  u  +  v  f(u) 

C  t 

vanishes  for  all  piecewise  smooth  closed  contours  C.  As  we  remarked  earlier. 
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one  need  only  verify  this  property  for  a  substantial  class  of  contours,  for 
example,  parallelograms  with  sides  parallel  to  two  fixed  directions. 

We  shall  begin  with  several  comments  which  will  facilitate  the  convergence 
proof. 

Conservative  Schemes.  Consider  a  conservative  scheme  based  on  a  centered 
diamond  shaped  stencil  with  a  generating  function  <p  derived  from  an  equation 
of  the  form 

H(u  ,u  )  -  H(u  ,u  )  +  G(u  ,u  )  -  G(u  ,u  )  =  0  , 
ne  ws  nw  es 

where  H  and  G  are  defined  by  (2.5)  using  normals  a  and  5.  Let  ft(i,j) 

denote  the  mesh  diamond  centered  at  (iAx,jAt)  and  regard  it  as  the  union  of 

four  congruent  triangular  regions  Tne*  Tnw»  Tws  abd  Tes  obtained  by 

intersection  with  the  four  standard  quadrants  centered  at  (iAx,jAt);  here  the 

triangles  share  a  common  vertex  and  have  hypothenuses  indicated  by  the 

subscripts.  Suppose  that  the  values  of  all  difference  approximates  under 

consideration  lie  in  a  small  neighborhood  N  of  a  fixed  state  u  e  Rn.  A 

simple  application  of  the  implicit  function  theorems  shows  that  if  the  matrices 

a  I  +  a  Vf(u)  and  0  1  +  8  Vf(u) 
t  x  t  x 

are  invertible,  i.e.  if  the  C-F-L  condition  holds,  then  there  exist  smooth 
maps  0  and  defined  in  a  neighborhood  of  (u,u)  such  that 

H(a,b)  =  at9(a,b)  +  axf{0(a,b)},  n(a,a)  =  a 
G(a,b)  =  0t^(a,b)  +  if»(a,a )  =  a  . 

Using  0  and  one  can  associate  with  a  grid  approximation  u(Ax)  of  the 
scheme,  an  everywhere  defined  piecewise  constant  function  u  =  u(x,t,Ax)  such 


if  C  is  any  closed  polygonal  arc  consisting  of  line  segments  joining  adjacent 

mesh  points:  if  we  put  u  *  0(u  ,u  )  in  T  ,  u  =  <|i(u  ,u  )  in  T  , 
r  ne  ne  n  w  nw 

u  *  0(u  ,u  )  in  T  and  u  =  il/(u  ,u  )  in  T  we  obtain 
w  s  ws  e  s  es 

/  v(u)ds  =  0  , 

3fl(i,j) 

for  all  mesh  diamonds  ft  and  consequently  (11.2)  by  addition.  Now,  if  u(Ax) 
satisfies  uniform  estimates  of  the  form  (10.2)  and  (10.3)  then  so  does 
u(x,t,dx)  and  it  follows  that  a  subsequence  converges  to  function  say  u  such 
that 

(11.3)  /  v(u) ds  =  0 

C 

for  all  parallelograms  C  in  t  >  0  with  normals  a  and  0.  Thus,  u  is  an 
exact  solution. 

Perturbations  of  Conservative  Schemes.  A  similar  argument  applies  to  the 
perturbed  schemes  introduced  in  Section  6  with  generating  functions  of  the  form 
+  q  where 

2  2 

| q(u  ,u  ,u  ) |  <  const. (|u  -  u  |  +  |u  -  u  |  )  . 

w  s  e  w  s  s  e 

Using  the  function  n  and  iji  associated  as  above  with  the  generator  <f>  we 
obtain 

(11.4)  /  v(u)ds  =  0(Ax)(!u  -  u  I'  +  |u  -  u  |2) 

3ft(i,j)  w  s  s  e 

for  the  corresponding  piecewise  constant  approximation  u  =  u(x,t,Ax)  obtained 

by  setting  u  =  0{u  ,u  )  in  T  ,  etc.  Therefore,  if  u(Ax)  satisfies 
n  e  ne 

uniform  estimates  of  the  form  (10.2)  and  (10.3)  it  follows  that  a  subsequence  of 
the  functions  u(x,t,Ax)  converges  to  a  function  u  satisfying  (11.3)  for  all 
associated  parallelograms  C  provided  that  the  maximum  magnitude  of 
interpolated  waves  vanishes  as  the  mesh  is  refined,  i.e. 
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lim  sup{|u(x,t)  -  u(x  +  Ax,t  +  At) |  +  |u(x,t)  -  u(x  +  Ax,t  -  At)|}  =  0 
Ax+O 

where  the  supremum  is  taken  over  all  mesh  points  (x,t),  say  in  a  given  bounded 
set.  Thus,  a  stable  and  convergent  sequence  of  difference  approximations  for  a 
perturbed  scheme  yields  an  exact  solution  provided  that  the  maximum  magnitude  of 
waves  in  any  bounded  set  of  t  >  0  vanishes  as  the  mesh  is  refined. 

The  Random  Choice  Scheme.  With  a  grid  function  u(Ax)  produced  by  the  random 
choice  method  we  shall  associate  the  standard  piecewise,  smooth  approximtion 
u(x,t,Ax)  as  follows.  Let 

R(a,b)  =  {(x,t)  :a-Ax<x<a+  Ax,  b  <  t  <  b  +  At}  . 

Let  u1  =  u.j(x,t)  denote  the  restriction  to 

(11.5)  R(iAx, jAt)  ft(i,j) 

of  the  solution  to  the  Riemann  problem  at  time  t  =  jAt  z  with  data 

(uw,un)  having  a  jump  at  x  =  iAx,  where  uw  and  uR  denote  the  values  of 
u(Ax)  at  the  west  and  north  vertices  of  Similarly,  let  u2  = 

u2(x,t)  denote  the  restriction  to 

(11.6)  R{ ( i  -  1 )Ax, ( j  -  1 ) At}  fl(i,j) 

of  the  solution  to  the  Riemann  problem  at  time  t  =  (j  -  1  )At  with  data 

(u.u  )  jumping  at  x  =  (i  -  1)Ax;  u  the  restriction  to 
w  s  3 

(11.7)  R{ ( i  +  1 )Ax,( j  -  1 ) At}  H(i, j) 

of  the  solution  of  the  Riemann  problem  with  data  (ug,ue)  at  t  *  (j  -  1)At, 
x  =  (i  +  1)Ax.  Define  u(x,t,Ax)  in  (2(i,j)  to  be  u1#  u2  or  u3  according 
to  the  location  of  (x,t)  in  one  of  the  three  corresponding  subregions  (11.5), 
(11.6)  and  (11.7)  (we  will  shortly  employ  the  definition  above  in  the  setting  of 
hybridized  schemes,  for  those  diamonds  Sl(i,j)  in  RC). 

Now  the  integral  around  the  boundary  of  a  typical  diamond  Q(i,j)  can  be 
expressed  in  terms  of  the  location  of  the  corresponding  sample  point  as  follows: 
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(11.8) 


where 


/  v(u)ds  =  E+(i  -  1 , j 
3ft(i,  j) 


-  1)  +  E  (i  +  1,j  -  1) 


E  =  f  u  -  R(x/At;u  ,u  )dx,  E  =  /  u  -  R(x/At;u  ,u  )dx 


Here  uw  and  ug  denote  the  values  of  the  corresponding  and  function  at  the 
west  and  east  vertices  of  ft(i,j).  Now  if  0  is  parallelogram  in  the  x  -  t 
plane  which  can  be  represented  as  a  union  of  mesh  diamonds  ft(i,j,Ax)  for  each 


Ax,  i.e. 


(11.9) 


ft  =u{ft(i,j,Ax)  :  (i,j)  e  A( Ax) } 


for  an  appropriate  index  set  A(Ax)  then 

l  E+(i  -  1,j  -  1)  +  E~(i  +  1,j  -  1) 

(i, j)eA(Ax) 

approaches  zero  as  the  mesh  is  refined  by  results  on  Liu  [26],  if  TVu0  k<  ^ • 

It  follows  that  the  limit  of  sequence  of  random  choice  approximations  converges 
to  an  exact  solution  u.  We  note  that  large  parallelograms  of  the  form  (11.9) 
exist  if  we  choose  for  example  a  sequence  of  mesh  lenghts  of  the  form 
(Ax)^  =  const. /2n.  For  general  mesh  lengths  a  trivial  modification  of  the 
above  argument  is  required  to  obtain  the  same  conclusion. 

Hybridized  Schemes.  Here  we  shall  associate  with  a  given  grid  function  u(Ax) 
a  piecewise  smooth  approximation  u(x,t;Ax)  by  following  the  procedure  above 
for  conservative  (or  what  is  the  same  perturbed  conservative)  schemes  if 
ft(i,j)  lies  in  C  and  the  procedure  for  the  random  choice  method  if  ft(i,j) 
lies  in  RC.  For  such  u(x,t,Ax)  the  contour  integral  around  the  boundary  of  a 
fixed  parallelogram  ft  of  the  form  (11.9)  can  be  expressed  as  the  sum  of  three 
terms  1^,  I2  and  I-j.  The  first  term  I1  records  the  contribution  from 
diamonds  ft(i,j)  in  C  due  to  the  perturbation  q: 

I1  =  0[Ax  £  {e2  :  e  crossing  a  diamond  ft ( i , j )  in  c}]  . 


L 
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The  second  term  I2  records  the  contribution  along  the  boundary  B  between  the 
two  regions  formed  with  mesh  diamonds  in  C  and  RC  respectively  and  takes  the 
form 

I2  *  0( Ax*TVu)  . 

B 

The  third  term  represents  the  random  choice  error  and  takes  the  form 
(11.10)  I3  =  I  {E+(i-1,j-1)  +  E~ (i+1 , j-1 )  :  ft(i,j)  £  RC  ft}  +  m(Ax), 
where  m(Ax)  represents  the  alterations  to  the  standard  random  choice  error  due 
to  the  use  of  the  modified  random  choice  generator;  as  we  showed  in  Section  7, 
m(Ax)  vanishes  as  the  mesh  is  refined. 

It  follows  immediately  from  (10.2)  that  I ^ ,  has  the  order  cr(Ax)  and 
thus  vanishes  in  the  limit.  We  shall  show  below  that  the  total  variation  of  the 

difference  approximations  u(x,t,Ax)  around  the  boundary  B  =  B(Ax)  has  the 

order  o(1/Ax)  and  thus  I2  vanishes  in  the  limit.  Finally,  the  argument  of 
Liu  for  the  random  choice  method  applies  with  an  only  trivial  modification  to 
show  that  I3  vanishes  in  the  limit. 

The  analysis  of  I2  proceeds  as  follows.  Fix  a  difference  approximation 
u  *  u(x,t,Ax)  and  a  time  strip  [0,T].  Let  K1  denote  the  set  of  all  mesh 
diamonds  in  C  [0,T]  such  that  either  (2  or  ft+  lies  in  RC  and  B^  the 

set  of  all  waves  e  with  the  following  three  properties:  e  is  incoming  with 
respect  to  a  diamond  (2  in  if  e  croases  the  ws  side  of  ft  then  ft" 

lies  in  RC;  if  e  crosses  the  se  side  of  ft  then  ft+  lies  in  RC.  Let 

K2  denote  the  set  of  all  mesh  diamonds  (2  in  RC  n  [0,T]  such  that  either 
(2  or  ft+  lies  in  e  and  let  B2  denote  the  set  of  all  waves  £  with  the 
following  three  properties:  e  is  an  incoming  wave  with  respect  to  a  diamond 
(2  in  Kjj  if  e  crosses  the  ws  side  of  (2  then  ft  lies  in  C,  if  e 
crosses  the  se  side  of  ft  then  ft+  lies  in  C.  We  have  B  =  B,,  u  B2.  For 
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simplicity  we  suppress  the  dependence  of  Bj  on  the  mesh  length.  We  shall  show 
that 

(11.10)  lim  AxTV  u(*,»,Ax)  =  0  , 

Ax*0  (Ax) 

a  similar  analysis  produces  the  same  result  for  B2>  To  this  end,  let  us  fix  a 
difference  approximation  u(x,t,Ax)  and  a  diamond  0  in  B^.  We  shall  analyze 
the  contribution  due  to  waves  crossing  the  ws  side  of  £2;  the  se  side  is 
treated  in  a  similar  fashion.  Let  us  therefore  assume  that  0  lies  in  RC 
and  for  concreteness  that  the  Riemann  problem  with  data  ( uw , ) ,  where  uw 
and  u@  denote  the  values  of  u  at  the  west  and  east  vertices  of  (2  ,  gives 
rise  to  two  shock  waves.  Let  and  (Y^»Y2>  denote  the  incoming  waves 

of  £2  crossing  its  ws  and  se  sides;  let  (a.^  ,a  )  and  <Jenote  the 

outgoing  waves  of  £2  crossing  its  wn  and  ne  sides.  We  consider  the 

following  cases. 

Case  1 .  Suppose  that  a  major  wave  terminates  in  £2  .  Then  the  total  strength 
S  of  waves  leaving  £2  and  entering  £2  is  less  than  ca.  If  we  associate 

with  £2  the  diamond  T  which  contains  the  initial  segment  of  the  major  wave 

-  2 

terminating  in  £2  we  have  H(T)  >  ca  and  therefore 

AxS  <  cAxo  <  cAxH(D/a 

Since  the  sum  of  squares  of  all  initial  segments  of  major  waves  is  finite,  i.e., 

I  H(D  <;  C 

it  follows  that  the  total  contribution  from  all  waves  of  the  above  type  (in  a 
fixed  strip  [0,TJ )  satisfies 

Ax  J  S  <  cAx  J  H(r)/o  <  cAx/a 
and  therefore  vanishes  by  (7.9). 
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Case  2.  Suppose  that  no  major  wave  terminates  in  ft  .  We  may  assume  that  ft 
contains  a  major  1-wave  and  that  a  1 -shock  crosses  the  wn  side  of  ft 

while  a  2-shock  crosses  the  ne  side  of  ft  .  Let  ft1  and  ft^  denote  the 
diamonds  whose  ne  and  wn  sides  coincide  respectively  with  the  ws  and 
se  sides  of  ft  .  We  shall  first  treat  the  subcase  where  ft1  and  ft2  both 
lies  in  RC.  Here  we  have 

e2  =  62  +  Y2  +  0D(n-)  • 

Now  if  y  lies  on  a  major  wave  then  ly^  >  a  and  we  obtain  an  estimate  of 
the  form 

|oP2l  <  1 Y152l  +  • Y, Y2 I  +  0{oD(ft“)} 

<  c{D(ft~)  +  . 

If  61  lies  on  a  major  wave  then  16^1  >  0  and  we  obtain 

1 a6 21  <  I «1«2I  <  clfft.,) 

But,  if  61  *  0  then  Y2  =  0  and  we  obtain 

I af?2 1  <  |o62|  +  0{oD(ft')}  <  c{D(ft“)  +  I(ft2)}  . 

Thus,  the  total  strength  of  waves  corresponding  to  the  subcase  where  ft1  and 
ft2  both  lie  in  RC  is  bounded  by  const. /o. 

The  subcase  where  ft^  and  ft2  both  lie  in  C  is  even  simpler  since  we 
have  an  estimate  of  the  form 


| B2 I  <  ca  <  cH(ft”)/o 
-  2 

using  the  fact  that  H(ft  )  >  ca  . 

If  tli  e  C  and  ft2  t  RC  we  proceed  as  follows.  If  Iy^  >  a  we  remark 

that 


|o3  |  <  c{D(ft")  +  I(ft3) 
as  in  the  first  subcase,  while  if  15^1  >  a  we  have 

I oB2 |  <  cH ( ft” ) 


as  in  the  second  subcase.  It  remains  only  to  treat  the  subcase  where 
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and  e  C.  Now  if  y^  lies  on  a  major  wave  then  the  analysis  of  the  first 

subcase  applies.  If  ^  lies  on  a  major  wave  then  we  proceed  as  follows.  Let 

p(Ax)  be  a  positive  function  such  that 

lim  p(Ax)  =  0  =*  lim  Ax/pc 
Ax+0 

We  distinguish  two  situatins  accordingly  as  10^1  *  and  |02I  <  per.  In  the 
first  situation  we  have 

Ax|02!  =  (Ax/|02I)02  <  cAx{52  +  y2  +  D(n“)}/po  . 

We  observe  that  1 6^ |  <  cK^l/a  using  the  fact  that 

15/21  <  clffl^  . 

But  we  need  only  treat  the  case  where  | 52 |  <  ca .  Therefore 

S2  <  clO^) 

and 

Ax|B2I  <  cAxtf^)  +  I*(fl")}/po 

since  y 2  is  associated  with  a  maximal  weight.  In  the  situation  where 
I B2 I  <  po  we  have 

ie2l  <  ulo^l 

and  consequently  the  sum  of  all  such  waves  B2  in  a  given  strip,  say 
[t,t  +  At]  satisfies 

I  Ib2i  <  ptvx 

where  X  is  the  associated  configuration.  Since  the  number  of  mesh  strips  in 
[0 ,T]  has  order  1/Ax  we  conclude  that  the  total  strength  of  all  such  waves 
02  in  [0,T]  satisfies 

I  1 62 I  <  cp/Ax  . 

This  completes  the  proof  of  (11.10). 

1  /2 

We  note  that  optimal  choice  of  p,  i.e.  p  =  (Ax/o)  ,  leads  to  a  growth 
estimate  of  the  form 
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1  /2 

TV  u(*,»,Ax)  <  const. /(Axa)  • 

B1UB2 

We  conclude  this  section  with  a  few  brief  remarks  concerning  the  first  term 
on  the  right  hand  side  of  (11.10).  In  this  connection  we  shall  first  comment  on 

the  argument  used  in  [26]  to  show  that  the  error  associated  with  the  random 

choice  method  vanishes  for  an  arbitrary  equidistributed  sequence.  Fix  e  >  0 
and  T  >  0.  In  [26]  the  time  strip  (0,T)  is  divided  into  substrips  of  equal 
length  k(e)  and  the  error  associated  with  each  substrip  is  estimated  using  a 
partitioning  of  elementary  waves,  cf.  Lemma  3.2,  [26].  It  is  then  shown  that  in 
the  limit  as  the  mesh  is  refined  the  sum  of  the  contributions  of  each  strip  has 
the  order  of  e.  The  partitioning  divides  the  set  of  elementary  waves  of  a 
given  random  choice  difference  approximation  into  two  classes.  The  total 
magnitude  of  waves  in  the  first  class  is  small  in  the  sense  that  their  total 

contribution  to  the  error  has  order  e  independently  of  whether  the  sequence  is 

equidistributed  or  not.  The  waves  of  the  second  class  in  a  given  substrip  can 
be  grouped  into  polygonal  arcs  Y  in  the  x-t  plane  whose  magnitude  and  speed 
of  propagation  are  "nearly  constant".  This  leads  to  an  individual  treatment  of 
the  arcs  Y  modelled  on  the  argument  used  for  a  single-wave  solution  to  a 
Riemann  problem.  Now  if  the  equidistributed  sequence  admits  a  rate  of 
approximation  in  the  sense  that  the  average  number  of  its  elements  with  index 
less  than  n  in  a  given  interval  J  C  (-1,1)  satisfies 

S(n,J)/n  *  |j|  +  c(J)<l>(n)  and  lim  ^ ( n )  =  0 

n+“» 

where  the  constant  c(J)  is  independent  of  n,  then  one  may  divide  (0,T) 
into  substrips  of  equal  length  k(Ax)  satisfying 

lim  k( Ax)  -  0 
Ax+0 
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and  again  add  the  contributions  to  each  stip  obtaining  in  the  limit  a  total 
contribution  on  the  order  of  e  provided  that  the  length  of  the  substrips 
approaches  zero  slower  than  the  reciprocal  of  the  error  ip,  i.e.  provided  that 

lim  \|;(  1/Ax)/k(Ax)  *  0  . 

Ax+0 

We  recall  that  the  hybridized  schemes  under  discussion  employ  an 
equidistributed  sequence  and  threshold  function  satisfying  (11.1).  Given  such 
i|»  and  o  one  can  employ  a  similar  two-class  partitioning  of  elementary  waves 
of  the  hybridized  shceme  to  strips  with  length  l  (Ax)  satisfying 

lim  <|>(1/Ax)/*(Ax)  =  0  and  lim  l(  Ax)/(aAx)1/2  =  0 
Ax+0  Ax+0 

and  show  that  the  total  contribution  to  the  error  from  all  waves  of  the  first 
class  plus  all  waves  of  the  second  class  which  form  arc  Y  with  length  greater 
than  say  i(Ax)/2  has  the  order  t  in  the  limit  as  the  mesh  is  refined.  Thus, 
one  need  only  show  that  those  segments  of  arcs  Y  with  length  less  than 
4(Ax)/2  passing  through  diamonds  in  RC  produce  a  total  contribution  which 
vanishes  in  the  limit.  To  this  end  let  us  fix  a  substrip  S(Ax)  of  (0,T) 
with  length  4(Ax)  and  consider  an  arc  Y  of  second  class  waves  e.  The  error 
associated  with  those  segments  of  Y  which  lie  in  the  set  ft  of  all  RC 
diamonds  in  S(Ax)  has  the  order  Ax|e|.  Without  loss  of  generality  we  may 
assume  that  the  magnitude  of  waves  of  Y  remains  a  constant,  denoted  |y|, 
since  from  the  results  of  [26]  the  deviation  from  a  constant  leads  to  terms 
whose  total  sum  vanishes  in  the  limit.  Hence  the  problem  is  to  estimate 

Ax  r(Y) | Y | 

where  r(Y)  denotes  the  number  of  segments  on  Y  which  correspond  to  diamonds 
in  £).  Now,  it  follows  from  our  previous  analysis  that  we  have  an  estimate  of 
the  form 
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I  T (Y ) | Y |  <  q{S(Ax)}/(cAx) 1/2  , 

Yes 

where  t(Y)  denotes  the  number  of  times  that  the  arc  Y  passes  from  ft  to 
c  c 

ft  or  from  ft  to  ft  and  where  the  sum  of  the  quantities  q  associated  with 
each  strip  S(Ax)  is  bounded  independently  of  Ax  for  fixed  T: 

£  {q{s( Ax) }  s  S(Ax)  C  (0 ,T) }  <  const  . 

Therefore,  the  contribution  associated  with  one  strip  S(Ax)  satisfies 

l  Axr(Y) | Y|  <  l  t(Y)*(Ax)|Y|  <  *(Ax)q{s(Ax)}/(oAx)1/2  . 

Yes  Yes 

Summing  over  all  strips  S(Ax)  in  (0,T)  yields 

£  £  Axr(Y)|Y|  <  const. £(Ax)/(oAx)1^2 

S  Yes 

which  implies  that  the  total  contribution  vanlahea  In  the  limit. 

12.  THE  ENTROPY  CONDITION. 

In  this  section  we  shall  show  that  the  limiting  solutions  u  of  our 
hybridized  schemes  satisfy  Lax's  entropy  condition  which  requires  that 

(12.1)  n(u).  +  q(u)  <  0  , 

t  x 

in  the  sense  of  distributions  where  n  is  a  strictly  convex  entropy  and  q  the 

associated  entropy-flux  [21].  Since  the  solutions  u  under  consideration  lie 
00 

in  BVn  1  it  is  sufficient  to  show  that 

(12.2)  /  v  n(u)  +  v  q(u)ds  <  0 

C  X 

for  a  substantial  class  of  contours  C,  e.g.  parallelograms  whose  sides  have 
slope  tAt/Ax.  As  a  corollary  of  the  entropy  condition  it  follows  that  the 
entire  sequence  of  associated  difference  approximations  u(Ax)  converges  to 
u  as  the  mesh  is  refined  in  those  circumstances  where  uniqueness  is  available 


-87- 


[9];  for  example.  If  the  initial  data  Uq(x)  gives  rise  to  a  piecewise 

Lipschitz  solution  u  to  a  genuinely  nonlinear  system  of  two  equations  in  the 

sense  of  [9]  and  if  u  satisfies  (12.1),  then  u  is  unique  within  the  class  of 
00 

BV  n  l  solutions  satisfying  (12.1).  The  class  of  piecewise  Lipschitz 
solutions  discussed  in  [9]  contains  in  particular  the  classical  solution  of  the 
Riemann  problem. 

Let  us  first  consider  an  arbitrary  conservative  scheme  in  class  K.  A 
simple  computation  using  the  compatibility  condition 

7n(u)7f(u)  m  7q(u) 

shows  that 


afcn{0(n,e)}  +  axq{8(n,e)}  -  atn{6(w,s)}  -  axg{9(w,s)} 

+  8tn{<|»(n,w) }  +  8xq{i|/(n,w) }  -  0tn{^(e,s)}  -  0xq{i)i(e,s) } 

=o(|u  -u|2+|u  ~  u  | 2 ) 

w  s  s  e* 

Thus,  the  everywhere  defined  piecewise  constant  approximations  u  *  u(Ax) 
associated  with  a  conservative  scheme  satisfy 


/  v(u)ds 

an 


0(Ax) ( |u  -  u  | 2  +  |u 
w  s  s 


uel2) 


for  a  typical  mesh  diamond  ft  with  values  uw,  ug,  ug  at  its  west,  south  and 
east  vertices.  Thus  a  convergent  sequence  of  difference  approximations 
u  =  u(Ax)  which  is  stable  in  the  total  variation  norm  yields  a  solution  u 
which  satisfies  (12.2)  for  all  parallelograms  C  with  sides  of  slope  ±At/Ax, 
provided  that  the  maximum  magnitude  of  waves  in  any  bounded  set  T  of  the 
x-t  plane  approaches  zero  as  the  mesh  is  refined,  i.e. 

lim  sup { | u( x  ±  Ax,t  ±  At)  -  u(x,t) |  :  (x,t)  t  T}  =  0 
Ax+0 


We  note  that  exactly  the  same  conclusion  can  be  drawn  for  any  scheme  whose 
generating  function  is  formed  by  a  quadratic  perturbation  of  a  conservative 


generator.  Of  course,  in  the  absence  of  a  condition  on  the  limiting  maximum 
wave  magnitude,  the  problem  of  determining  whether  or  not  the  entropy  condition 
is  satisfied  is  substantially  more  difficult. 

For  the  original  version  of  the  random  choice  method,  the  entropy  condition 
was  verified  by  Lax  [21].  For  the  deterministic  verion  involving  an 
equidistributed  sequence  it  follows  immediately  from  the  results  of  Liu  [26] 
that  the  entropy  condition  is  satisfied.  In  the  latter  case  one  has 

/  v(u)ds  <  E+(i  -  1 , j  -  1)  +  E“(i  +  1,j  -  1) 

3fl(i, j) 

where 

Ax  0 

E  m  [  n(u  )  -  n(R(x/At;u  ,u  ))dx;  E  *  /  n(u  )  -  n(R(x/At;u  u  ))dx 

'  n  we  '  n  we 

0  -Ax 

where  uw  and  ufi  denote  the  values  of  the  corresponding  grid  function  at  the 
west  and  south  vertices  of  the  mesh  diamond  The  inequality  in  (12.3) 

arises  from  the  fact  that  the  discontinuities  of  u(Ax)  in  fl(i,j)  satisfy  the 
entropy  condition  and  the  sum  of  all  terms  of  the  form  E*  vanish  in  the  limit. 

Therefore  in  view  of  our  previous  analysis,  it  follows  that  the  limits  of 
our  hybridized  scheme  satisfy  the  entropy  condition.  We  need  only  remark  that 
the  results  of  Section  11  show  that  the  interior  boundary  contribution 
(associated  with  the  set  C  and  RC)  to  the  contour  integral 

/  v  n(u)  +  v  q(u)ds  , 

c  r  x 

vanishes  in  the  limit,  since  the  quantity  Ax  times  the  total  variation  of 
u(Ax)  around  said  boundary  vanishes  in  the  limit. 
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